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ABSTRACT 


This thesis details the system identification and initial system validation of the an 
Ultra-Quiet Vibration Isolation Platform (UQP). With the move toward lighter and more 
flexible spacecraft, the effects of vibration are of immense concern. As natural or passive 
damping becomes less effective in controlling undesired vibrations, active vibration 
control becomes essential. The UQP uses a special configuration of the six degree of 
freedom Stewart Platform with piezoceramic strut actuators and geophone sensors. This 
combination gives an extremely sensitive and responsive six degree-of-freedom active 
vibration control system. Each actuator was designed to be controlled independently 
without coupling with other actuators. In order to develop control laws, the plant must be 
identified in terms of system zeros and poles and the uncoupled design validated. 
Dynamic modeling using parametric estimation methods can accurately describe a 
complex system. Using parameter estimation methods, models of the actuator system 
dynamics were obtained. A simple lead-lag controller was applied to individual actuators 
then all six actuators acting simultaneously to verify system coupling. Significant 


interaction between base adjoining actuators was discovered. 
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I. INTRODUCTION 


A. PURPOSE FOR RESEARCH 


With increased use of lightweight materials, spacecraft will become more flexible 
and susceptible to excitation from sources of vibration. With loss of stiffness and rigidity, 
passive damping alone may be insufficient to attenuate disturbances to an acceptable 
level. Active Vibration Control (AVC) is rapdily emerging as a method of eliminating or 
reducing unwanted vibration. The purpose of this thesis is to provide a detailed 
description and system identification of the Ultra Quiet Vibration Isolation Platform 
(UQP). It will be the basis for future research in the area of AVC. The UQP is an 
excellent sytem for AVC and can be applied to a broad range of disturbances. A 
disadvantage of the UQP is the complexity of the plant dynamics and kinematics. 
Accordingly the complexity of the control system design is increased. This thesis will 


charaterize the plant for future control design on the UQP. 


B. SCOPE 


The “plant” is the UQP which consists of six piezoceramic control actuators. This 
thesis will use parameter estimation methods to create a dynamic model of the UQP 
contro] actuators with the ultimate goal of extracting actuator zeros, poles and transfer 
functions. These actuators are designed to operate independently. Theoretically, each 
actuator can be controlled without coupling or interaction with the other actuators. This 
assumption is of extreme importance as it critically impacts upon the performance of the 
system and control system design. This assumption will be checked for validity by 
applying a lead-lag feedback compensator to an actuator and then to all six actuators 


simultaneously. 





Il, EXPERIMENTAL SETUP 


A. BACKGROUND 


In many spaceborne applications the dynamics and control of vibrations must be 
addressed as a multiple degree-of-freedom (DOF) problem. Translations and rotations 
about all three axis must be considered. The Stewart Platform is the ideal mechanism for 


multiple DOF vibration control applications. 


1. The Stewart Platform 


The original motivation put forth by Stewart {Ref. 1] was to design a 
mechanism capable of simulating flight conditions to train pilots. In order to be realistic 
it had to be capable of translating and rotating in three directions just as a real aircraft. 
The original configuration consisted of a triangular plate and a rigid parallel base 
connected by six legs in a “linear coordinate” leg system. At each triangle vertex two 
legs were attached in a three DOF joint mechanism. These legs, with hydraulic actuators, 
were mounted to two DOF joints at the base. The advantages of choosing this 
configuration were inherent rigidity and absence of bending moments. Additionally with 
this configuration the only forces present were in the plane of the leg. A similiar six leg 
arrangement had also been used by a machine designed to study tire to ground forces [Ref 
2]. In this system computer controlled jacks were used as actuators. A cubic 
arrangement was devised such that the relationship of one actuator to any other is the 
same (see Figure 2.1). This arrangement showed inherent stability and the capability for 


accepting large moments. 





Figure 2.1 
Cubic Arrangement 


Stewart Platform mechanisms have been the subject of studies as multiple DOF parallel 
systems [Ref. 3 and Ref. 4:p 46]. The Stewart Platform configuration has shown wide 
applicability from motion simulators to robotics and now Active Vibration Control. The 
UQP employs the cubic configuration of the Stewart Platform which has an extremely 
important property with respect to AVC. With the exception of inertial loads and gravity 
forces, all other forces are carried axially. The significance is that if axial forces due to 


2 


vibration can be eliminated the vibration is eliminated [Ref. 4:p. 46]. 


B. HARDWARE CONFIGURATION 


The experimental hardware is divided into two sections: the UQP and the Power, 


Control and Analysis Hardware. 


1. Ultra-Quiet Vibration Isolation Platform (UQP) 


The UQP is an adaptation of the Stewart Platform for AVC designed and 
built by CSA Engineering INC. It consists of a circular plate attached to a rigid base by 


six variable length actuators (see Figure 2.2). 
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igure 2-2 


Ultra-Quiet Vibration Isolation Platform 


The base consists of an aluminum plate supported by aluminum stringers and longerons. 
The base is used to simulate a rigid spacecraft. Mounted inside one of the “spacecraft” 
bays is an AURA Bass Shaker which is used as a known disturbance source. The entire 


aparatus is bolted to a 38001b NEWPORT optical table. 


a) Control Actuator Struts 


The UQP uses six mutually orthoganol struts to provide control 
over six degrees of freedom. Each strut consists of a piezoceramic actuator and a 


geophone sensor. The use of piezoceramics for shape and vibration control is becoming 


increasingly widespread. For interested readers, further information on piezoceramics 
and their applications in vibration control can be found in Reference 8. By taking 
advantage of the cubic configuration, all six struts can be considered linear motion 
actuators [Ref 1]. Stewart Platform based AVC mechanisms using the same cubic 
configuration with Terfanol-D actuators [Ref. 4] and voice coil actuators [Ref. 5] have 
also been developed. These platforms displayed significant reduction of vibration over 
passive means alone. Figure 2.3 is a basic diagram of each actuator. The actuator stroke 
is approximately 50 microns. In the passive (open loop) mode it provides moderate 
damping to low frequency vibrations. In the active (closed loop) mode it provides 
damping to higher frequency vibrations. Passive damping is provided by a flexure with 
damping material. An extremly sensitive geophone measures velocity. This velocity is 


the error signal fed into the control Jaw. 








Figure 2.3 
Actuator 


In general applications, geophones consist of wire coils supported by soft springs under 
the influence of a magnetic field. Vibrations cause movement of the magnet relative to 
the stationary coil inducing a voltage proportional to velocity [Ref. 6:p 449]. The 
stiffness K refers to the stiffness of the passive stage. The stiffness K, corresponds to the 


piezoceramic stack actuator (PZT). 


b) Disturbance Shaker 


To create a disturbance source against which the performance of 
the UQP and applied control can be measured, a known disturbance source was mounted 
to the base. A sinusoidal waveform from a signal generator is amplified and input to a 
AURA Bass Shaker (see Figure 2.4). This simulates a cyclic disturbance. It has a 
resonant frequency at 42 Hz which is used as the disturbance frequency in initial 


experiments. 


Disturbance 
Shaker 





Figure 2.4 
Disturbance Shaker 


Z: Power, Control and Analysis 


To operate the UQP and extract meaningful experimental data requires 


several important components (see Figures 2.5 & 2.6). 


Control] Input-Sensor Output 


Sensor Output- Channel of Interest 


Ld 
XX 
aaar-<) 





Disturbance Shaker 


KEPCO Operational 
Amp. & Pwr. Supply 










HP 33120A Function/ 
Arb. Waveform Generator dSPACE 





Figure 2.5. 
Power, Control and Analysis Hardware Configuration 





Figure 2.6 
High Voltage Power Supply (HVPS)(top) and Active Vibration Control System (AVCS) 
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The current configuration requires a MATLAB capable PC with an RS232 Y/O port. The 
control design is done in MATLAB. Once the design is completed, the MATLAB code is 
translated into C code. The C code is converted to machine language and down loaded 
via the I/O port to the CSA Active Vibration Control System (AVCS) (see Figure 2.6). It 
provides the interface for control implementation and geophone sensor output. Its main 
component is a digital signal processor (DSP). Once the machine code has been loaded 
the AVCS applies the control to the actuators. It receives the feedback signal from the 
geophones, processes and feeds it into the control algorithm. The output is sent to a HP 
35665A Dynamic Signal Analyzer for analysis. The CSA High Voltage Power Supply 
(Figure 2.6) supports the power equirements of the piezoceramic actuators and is not 


required for use in the open loop mode. 
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lil. SYSTEM IDENTIFICATION 


In order to design an effective compensator or controller for a dynamic system the 
plant must be well known or characterized. By the use of system identification 
techniques a plant model can be obtained. The use of a model is extremely beneficial in 
cases when the dynamics of a plant are not well known or complex. Although articles 
have been written investigating the dynamic and kinematic complexities of Stewart type 
platforms [Ref. 3.], active vibration control using Stewart Platform based mechanisms is 
a relatively new field. Experimental dynamic modeling can be used to bypass the need 
for a complete theoretic dynamic analysis of the UQP. The goal of identifying the UQP 
plant was to obtain an accurate plant transfer function which can assist control design. Of 
particular interest were system modes or resonances and interaction or couplings among 
actuators. An initial key assumption on the actuators is that they each act independently 
(uncoupled) and can be controlled as such. Validation of this assumption is one of the 


goals of identification of each plant. 
. DYNAMIC MODELING 


There are two categories of system identification used in dynamic modeling: Non- 
parametric and Parametric. Non-parametric methods are used to estimate transient 
response, spectra and frequency functions in order to gain basic knowledge of the system. 
Parametric methods are used to obtain the mathematical parameters or elements used by 
well know system modeling algorithms. These models are broken down into two 
subcatagories. “Tailor Made” models are built based on physical principles (i.e. the 
parameters have some type of physical interpretation). “Ready Made” or “Black Box” 
models are general in nature and use parameters that have no direct physical 
interpretation. They are used to characterize the relationship of output to input of a 
system. By treating the UQP as a “Black Box” and analyzing the input-output 
relationship it is possible to obtain an accurate model of the system by using a Ready 


Made model. [Ref. 10] 


1. Background 


The concept behind the use of transfer functions is to describe the 
relationship between the input and output of a system and can be expressed in terms of 
either continuous or discrete time. The use of discrete time representations are preferable 


when the data is sampled in nature. 


unction 


Figure 3.1 
Input-Output Relationship 





One of the most important tools used to gain an understanding of the behavior of a system 
is by looking at its impulse response. For a linear system (plant) the impulse response 1s 


the output when the applied input 1s an impulse at t=0. 


O(t) y(t)=h(t) 


Figure 3.2 
Impulse Response Relationship 
If the impulse response h(t) of a system is known, the output to a given input can readily 
be determined by convolution of the input with the impulse response: 
y(t) = u(t)* h(t) (3.1) 


0) 


y(t)= Fh(t)eu(t—m) (3.2) 


m=-—oo 


The impulse response is a key element of non-parametric system identification. 

An essential element of system identification is characterizing the output. 
In an ideal system the output would simply be the input modified by the system transfer 
function. In physical systems however, the output also contains elements resulting from 


internal and external disturbances (see Figure 3.3). 
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Disturbance Disturbance 


Plant=h(t) 


y(t)=Total Output 
w(t)=Undisturbed Output 
o(t)=Disturbance term 





u(t) y(D=y(H+o(0 





Where: 


Figure 3.3 
In order to develop a model of the input output relationship, define: 
y(t)= y(t) + o(¢) 
y(t) = G(q, p)u(t) (B3.achic) 
o(t)= H(q, p)e(t) 


where €(t) is white noise, p is the parameter vector: 


Py 
p=|: (3.4) 
P,, 
and qr is the shift operator based on sampling interval 7: 
97 V(t) = y(t + T) 
Gz y(t) = y(t-T) 
Introduce the general difference equation: 
Magee) )-a,yi— 27) peta ye NT )— bu) + but —17)+...+b,u(r— LT) 
(3.6) 


(3.5 a,b) 


where T is the sampling interval of the discrete time representation. This standard 
expression 1s a specialized geometric series that relates input and output of a system using 
past values of input, output and a set of parameters a; and bj. By applying the shift 


operator qr, equation (3.6) becomes: 


= ig) ze a ae - 
y(|+a,97" +a2qT +° +a,97" | = u(0)| bo +bgyp +boqz + +b,g7" | 


= = = 
y(t) | bo thar +6247 +b ng” | 


1 lI+aiqp tangy" gre +ay97" | “_ 
Define: 
Ga,p)= 2" 3.8) 
and introduce a time delay of nk samples on the input: 
Que = [boar thay) thogp +t dpa | (3.9) 
A(q) ll+a,97" age +Onq97" | 


Changing variables to be consistent with the notation used by Ljung and Glad [ref. 11] 
and the MATLAB System Identification Toolbox (SITB) (Ref. 12]: 


Gia.) = B® boar +hyg*) +bygh >) +. -4b, ge) | ait 
)p)= = end 
F(q) I+ faz" +foqr tot fap gp 


Applying to equation (3.3b) and equation (3.6): 
WD) ate ieee) +S nf yw(t—nfT) = bou(t —nkT)+---+by,pudt—(nb+nk)T (3.11) 


Similarly the disturbance term becomes: 


S] =9 a 
C(q): lI+c197 +C24T +tencG7™ | 
HG, )) == 


a ee OS Co 
D(q) ll+d\gp!+dqqp + -+dng a7” | 


The coefficients b;, c;, dj and f; comprise the parameter vector p. These .coefficients 
represent the unknown system parameters. Essentially they approximate a complex 
system’s physical parameters but have no direct correlation to any specific physical 
quantity. The parameters nb, nc, nd and nf characterize the order and type of ready made 
model. The model constructed by equations (3.3) through (3.12) 1s shown in equation 


(3.13). It represents the Box-Jenkins model! [Ref. 10]. 
y(t) = G(q, p)u(t) + A(q, p)et) (Colley 


' Named after the statisticians G.E.P. Box and G. M. Jenkins 
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The following models are variations of the Box-Jenkins model. The most 
simple one is the case where the disturbance is not modeled. 
y(t) = G(q, p)u(t) + E(t) (3.14) 
Equation (3.14) 1s known as the Output Error method. The difference between the actual 
output and undisturbed output is manifested in the white noise term e(t). The next model 
utilizes the same poles for the input and disturbance (noise) models defined as: 
F(q) = D(q) = A(q) = 1+ agp +++ Gna Gp” (3.15) 
Applied to equations (3.10), (3.12) and (3.13) yields: 


TOE Ay ucty+ We (3.16) 
OF: 
A(q)y(t) = B(q)u(t)+ C(qg)e(t) (erly) 


This 1s known as the ARMAX (AutoRegression Moving Average eXtra input) model. 
The final model that will be described is the AutoRegression eXtra input(ARX) model 
(equation 3.18). [Ref. 10] 
A(q) y(t) = B(q)u(t) + E(t) (3.18) 
The next step is to describe how these ready made techniques model a 
system based on past values of output and input. Essentially these algorithms attempt to 
predict the output based on a given set of input-output data. Introduce the concept of 
One-Step-Ahead prediction of output. From equations (3.6) and (3.14), we have: 
y(t) = —a, y(t — 1)-...-a_, y(t — na) + b,u(t — nk)+...+b,,u(t—nk —nb+1)+e(t) (3.19) 
Since €(t) is white noise and cannot be predicted, the prediction becomes: 
3 (t, p) = —a,y(t-1)-...-a,, y(t —na) + b,u(t—nk)+...+b,,u(t—nk—nb+1) (3.20) 


Expanding to the general case of equation (3.13) and dividing out the noise transfer 


function: 
H™ (q4,p)y(t) = H~ (4, p)G(q, p)u(t) + e(t) (3.21) 
y(t)— y(t) + H"' (4, p) y(t) = Hq, p)G(q, p)u(t) + E(t) (3.22) 
y(t) +[-1+ H™(q, p)]y(t) = A” (4, p)G(q, p)u(t) + E(t) (rz2) 
y(t)=[1- H"(q, p)ly(t) + A (q, p)G(q, p)u(t) + E(t) (3.24) 


the prediction becomes: 
¥(t,p)=[1-H"(q, p)ly() + H'(4, p)G(q, p)u(t) (3.25) 
Since equation (3.25) only contains past values of the output y(t) and input u(t) the 
difference between the prediction and the output (the prediction error) is: 
E(t, p) = y(t)— y(t, p) = E(t) (3.26) 
where €(t) 1s white noise. [Ref. 9: p. 56 and Ref. 10: p.235] 


B. FREQUENCY RESPONSE OF ACTUATORS 


The diagram in Figure 3.4 is a modification of the experimental set up in Figure 


2.5 used to extract the input-output data. 
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Figure 3.4 
Experimental Setup 
The purpose behind measuring the frequency response of each actuator 1s to obtain both 
basic insight into the plant and a model verification baseline. Each actuator was excited 
by a50 mV “pink noise” source provided by the HP Dynamic Signal Analyzer (HPDSA). 
The “pink noise” is random noise which provides 3dB roll off per octave. The motivation 
to use pink noise 1s to place an equal amount of energy in each octave band. The 


response is sensed by the geophone of the excited channel (actuator) and fed into the 


HPDSA. The HPDSA computes and plots the time averaged frequency response which 
are shown in Figures 3.5 through 3.10. 
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Figure 3.5 
Frequency Response of Actuator #1 
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Frequency Response of Actuator #2 
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Frequency Response of Actuator #3 
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Figure 3.8 
Frequency Response of Actuator #4 
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Figure 3.9 
Frequency Response of Actuator #5 
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Frequency Response of Actuator #6 
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The plots show that the plant transfer functions for the actuators are basically the same. 
The resonant dynamics of the struts are readily discernible, specifically the actuator 


natural frequency at 1.4 kHz. 


C. MODEL STRUCTURE SELECTION 


The type of system and method of data collection are major factors determining a 
suitable model type for a given application. For example, the Output Error, OE, would be 
best applied when the disturbances (which are not modeled) are small. The ARX is the 
most widely used because it is the simplest and serves as a good baseline to select a 
model. However, it uses the same poles to model noise as those used to model the 
system. This can create problems in system identification requiring the use of higher 
order models than would otherwise be necessary. The Box-Jenkins method models the 
system dynamics and disturbances separately with no common parameters. The ARMAX 
also models the disturbance but uses the same poles as the system dynamic model. It 
differs from the ARX by the addition of the C(q) polynomial. The Box-Jenkins and 
ARMAX both provide detailed system models and can be extremely accurate. The 
decision on which to use is based on the nature of the disturbance. If “noise” enters the 
process in the early stages or is carried through the plant (see Figure 3.3) the ARMAX is 
preferable over the Box-Jenkins model which is more adept to characterizing noises that 
enter later in the system. In the case of the UQP, the ARMAX would appear to be the 
preferred method if there is interaction between actuators because the resulting 
disturbances would enter into the process early on. Additionally the disturbance caused by 
another actuator should show similar dynamics as the plant under experimentation. The 
flow chart in Figure 3.11 represents the process used in building a model of each actuator. 
As Figure 3.11 shows, the process of system identification and modeling is a step-by-step 
iterative process. The basic idea was to develop and validate a basic ARX model structure 
(i.e. na, nb and nk) for all six actuators and then apply it to the more complex ARMAX 
and Box-Jenkins methods and determine the best model. The application of this process 


on actuator number one will be detailed in the following sections. [Ref. 10] 
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Figure 3.11 
Model Selection Flow Chart 
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The computer program used to accomplish this can be found in appendix A. 


1. Data Collection 


The process of system identification and dynamic modeling relies on a set 
of plant input-output data. For the UQP each actuator was provided an excitation and the 
response of all six actuators (as measured by each respective geophone) recorded (see 
Figure 3.4). Using the HPDSA as the excitation source, 50mVRMS “Pink Noise” was 
applied to the actuator via the 25 pin connection on the front of the AVCS (see Figures 
2.6 and 3.3). Under normal operation the ribbon cable on the front of the AVCS is 
connected to close the control loop. By disconnecting the cable, each actuator can be 
accessed individually. The geophone sensor output was connected via a BNC interface 
box to a dSPACE system. The dSPACE was used to gather, display and save all six 
channels of data simultaneously. For comparison purposes the active channel sensor 
output (actuator under excitation) and excitation source was also input to the HPDSA. 


For each actuator the output corresponding to the excitation source was measured at a 


2) 


sampling frequency of 10 kHz for two seconds. This gave 20,000 input/output sample 
points. The choice of sampling frequency was based on the desire to obtain a sufficient 
number of sample points from which to construct and validate a model. The first 10,000 
(input and output) sample points are used by the parameter estimation algorithms to build 
the model. The last 10,000 points (of input) are used in model validation. If the sample 
points used to build the model are also used in validation, the true accuracy and validity 


of the model will be corrupted. 


a) Data Preparation 


Prior to using the data to build a model it 1s necessary to pre-treat 
the data. The most common factors adversely affecting collected data are: 
1. Drift, offset, trends and low frequency and/or periodic disturbances. 
2. Outliers or faulty data points. 
3. High frequency disturbances. 
Before treating for these possible problems the data is given a preliminary analysis to see 
if they do in fact exist or could be the possible source of erroneous results [Ref. 9:p 386]. 
External sources are the main cause of drifts, trends and low 
frequency disturbances that are either impossible or undesirable to model. This can be 
remedied by removing external disturbances, using the noise model to account for the 
disturbances and/or using an algorithm to de-trend the data. Due to the extreme 
sensitivity of the geophone sensors it is impossible to remove all the external 
disturbances. The latter two options were selected. It is assumed that the external 
disturbances received by the geophones are relatively small compared to those entering 
into the process early on and will be adequately addressed by the noise model. The offset 
problem is a result of the failure of input-output data to correspond in a consistent 
manner. This causes the model to waste parameters in attempt to adjust these levels. 
This is compensated by a MATLAB function included in the computer code (see 
appendix A.). [Ref. 12:p1-63] 


Le 


In any data collection effort there are obviously erroneous values. 
However, removal of these “outliers” requires caution. Removal can either be manual 
(i.e. selecting bad values by hand) or by some type of automatic algorithm. The best way 
to determine if outliers are problematic is to check the model residuals for excessively 
large values that stand out. This method was chosen due to the vast amounts of sampled 
data used. [Ref. 12:p1-63] 

The final source of data deficiency is the presence of high 
frequency disturbances outside the region of interest. In this case the primary concern is 
the performance of the actuators up to and including the resonance frequency. A 10th 
order Butterworth filter with a pass band below approximately 1.5kHz was used to 


eliminate disturbances outside the range of interest [Ref. 11: p 272]. 


2: Delay and Parameter Number/Structure Selection 


Prior to applying one of the aforementioned parametric models to the 
gathered data, the input delay nk, and the structural elements na, nb, nc, nd and nf must 
be determined. This was done using the Model Structure Selection functions in the 
MATLAB SITB. The SITB has functions to calculate ARX-based and OE-based 
structures. Figure 3.12 represents a flow chart of the algorithm used to select the delay 
and structure. Since the assumption was that the ARMAX model would provide the best 
model for the UQP, an ARX-based structure selection algorithm was applied (see 


appendix A). 
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Structure Selection Algorithm 


The process of selecting na and nb as well as the number of total parameters from which 
they are derived is an iterative process. Since the input delay should be the same no 
matter what the order of the model is, a second order ARX model will serve as the 
foundation of this process. Table 3.1 shows the results of nk values ranging from 1 to 10 


applied to a second order ARX model. 
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Table 3.1 
Input Delays 








The numbers represent the logarithms of a quadratic loss function based on the selected 
nk for na=nb=2 (second order model). Once an input delay has been selected a plot of 
the number of parameters used versus loss function is produced. Based on the results in 
Table 3.1, the best choice for an input delay for actuator number one appears to be nk=2. 


Figure 3.13. is the resulting plot. 
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Figure 3.13 
Parameters Vs. Loss Function for Actuator #1 nk=2 


The multiple results for a given number of parameters is due to the fact that there can be 


several different structures (1.e. combinations of na and nb). The SITB has functions that 


os 


automatically compute the structures based on the minimization of Akaike’s Information 
Theoretic Criterion (AIC) and Rissanen’s Minimum Description Length (MDL)[Ref. 11]. 
After a structure has been selected it can be applied to a parametric estimation method to 
develop a model. Even with the availability of functions capable of selecting the 
optimum number of parameters based on minimization techniques, it can be seen in 
Figure 3.13 that there is a point of diminishing returns beyond which the addition of more 
parameters is of little benefit. Further, the addition of more parameters than necessary 


will actually cause a loss in model accuracy. Referring to Figure 3.13 the structures 


chosen are listed in table 3.2. 


Number of Parameters Structure [na nb nk| 


Table 3.2 
Structures For Actuator #1 Input Delay nk=2 





The notes column provides information on how the structure was obtained. AIC and 
MDL based structures were automatically generated along with the plot in Figure 3.13. 
The term “Automatic” means the number of parameters in the first column were manually 
extracted from the plot in Figure 3.13 and entered in to an SITB function that 
“Automatically” selects the structure. The term “Manual” refers to the complete structure 


being manually selected. 


3: Model Structure Application and Evaluation 


After determining the ARX model structure for actuator number one, a 
model was computed using the algorithm in appendix A. There are numerous methods 
for measuring the accuracy and validity of a model. The ones utilized to validate the 
actuator models are described as they were applied to the ARX model of actuator number 


one using the structure [na=15 nb=15 nk=2]. 
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models are described as they were applied to the ARX model of actuator number one 


using the structure [na=15 nb=15 nk=2]. 


a) Frequency Response 


One of the quickest ways to determine a model’s validity is to 


compare the Bode plot of the system to that of the model. 
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Figure 3.14 
Frequency Response of ARX[15 15 2] Model (Actuator #1) 
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The model’s frequency response is similar to that of the actual system (see Figure 3.14). 
The areas of critical concern are the model behavior at resonance and zero dB crossing for 
magnitude and -180 degree crossing for phase. The phase is reasonably close however, 
the magnitude is significantly different. The general shape and resonance peak are in 


agreement but there is a distinct 10-20dB difference between the actuator and the model. 


b) Output Comparison 


In evaluating a model, as stated earlier, it is important to use 
different input sample data to insure an accurate assessment of the model. In modeling the 
actuators a 10,000 sample model validation set was held aside for this puppies. Applying 
the validation input to both the model and the system and comparing the outputs, gives an 
accurate measure of the model’s validity and how close it truly simulates the actual 
system. Figure 3.15 compares the model output to the actuator output in response to the 


same input. 
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Cc) Auto and Cross Correlation Functions 


Figure 3.16 provides plots of the auto correlation function of 
residuals and cross-correlation function between the input and residuals from the output. 
The difference between the model output and actual output is called the residual. 
Basically, it is what is left unaccounted for by the model. Recalling equation (3.26): 

E(t, p) = y(t)— y(t, p) = e(t) 
where €(t) represents the residuals. If €(t) is purely “white noise” it is independent of the 
input. If there is correlation between e(t) and u(t), there are elements in the output from 


the input that are not explained by the model. 
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Figure 3.16 
Auto and Cross-Correlation Functions for ARX [15 15 2] Model for Actuator #1 


The horizontal lines in Figure 3.16 indicate a 99% confidence level. If the function 
crosses these lines, there is a correlation between e€(t) and u(t) at that point. From Figure 
3.16 there appears to be correlation between e(t) and u(t-/). This indicates that the input 
delay might need to be modified to one. The auto-correlation function is acceptable. A 


plot of residuals versus sample number for 100 samples is given in Figure 3.17. 
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Residual of ARX Model [15 15 2] 


Residual 





-0.1 
5000 5020 5040 5060 5080 5100 
Sample 


Figure 3.17 
Residuals Vs. Sample ARX[15 15 2] For Actuator #1 
Over this range of samples the residuals are relatively low. For an average output 
magnitude of approximately 0.18, the average residual magnitude is on the order of 0.03 
or roughly 17% of the output magnitude. Although this percentage is a very rough way of 
comparison and cannot be considered by itself to determine model validity, it does give a 


basic idea of how much is left unexplained by the model. A much better method is to use 


the mean square fit. [Ref.10] 


d) Zero-Pole Plot 


The final evaluation tool is the zero-pole plot (see Figure 3.18). 
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Figure 3.18 
Zero-Pole Plot of ARX [15 15 2] For Actuator #1 


The ellipses correspond to confidence limits to three standard deviations [Ref.12]. This 
plot shows that there are numerous possible zero-pole cancellations indicating that a 


model of lesser degree is desirable. 


4. Model Application and Evaluation (Data Prefiltered) 


Based on analysis the ARX [15 15 2] is not a very good model. Following 
the flow chart the same analysis will be repeated on data pretreated with a 10" order 


Butterworth filter. 
a) Frequency Response 


By prefiltering the data the correlation between the model 
frequency response and the actual output frequency response is improved. However, 


there is still an approximate 10 dB difference between the two responses. 
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Figure 3.19 
Frequency Response of ARX[15 15 2] Using Prefiltered Data (Actuator#1) 
b) Output Comparison 


Figure 3.20 shows an improvement in the correlation between the 


model and actual output. 
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Figure 3.20 
Prefiltered ARX[15 15 2] Model Output Vs. Actual Output 


Cc) Auto Correlation and Cross Correlation Functions 


As a result of prefiltering there is a degradation of the auto and 
cross correlation functions (see Figure 3.21). However, the auto-correlation of residuals 
is relatively small. For cross correlation functions where the crossing occurs in negative 


lag, this is an indication that the data was possibly collected during feed back vice a 


problem with the model [Ref. 10:p 284] 
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Figure 3.21 


Auto and Cross Correlation Functions of Prefiltered ARX[15 15 2] (Actuator #1) 
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d) Residuals 


Comparing Figure 3.22 with Figure 3.20, the former residual plot 
is misleading. This anomaly occurs when large order models (large number of 
parameters) using prefiltered data are applied to the SITB function. This function uses 
statistical means to calculate the prediction error and the residuals. A simple plot of the 


actual output minus the model output is shown in Figure 3.23. 
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Figure 3.22 
Residuals of Prefiltered ARX[15 15 2] (Actuator #1) 
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Pigure 3.23 
Actual Output Minus ARX[15 15 2] Output (Actuator #1) 
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The average “residual” magnitude for untreated and pretreated data are 0.0041 and 
0.0032 respectively, or 2% and 1.7% of the average output magnitude. As previously 
stated, these are very rough measures of the accuracy but there is a clear improvement in 


the model by prefiltering the data. 


e) Zero-Pole Plots 


The zero-pole plot of the prefiltered model also shows 


improvement over the unfiltered model (see Figures 3. 18 and 3.24). 


Zero Pole Plot ARX Model {15 15 2] Model *Data Prefiltered* 





Figure 3.24 
Zero-Pole Plot For Prefiltered ARX [15 15 2] (Actuator #1) 


However, there are still numerous possible zero-pole cancellations, indicating that use of a 


lower order model is necessary. 


a: Model Structure Acceptance for Actuator Number One 


The goal of the model structure application and evaluation process is to use 
the ARX based technique to find an acceptable model structure that can be applied to 
more accurate ARMAX and Box-Jenkins methods. Knowing that prefiltering data 


provides a better model, the process is repeated. Through the first iteration it was also 
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determined that a delay change would possibly be beneficial. Retuming to the delay 


selection process, Figure 3.25 plots the loss function for a given number of parameters. 
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Figure 3.25 


Number of Parameters Vs. Loss Function For nk=1 (Actuator #1) 


The structures evaluated are contained in table 3.3. 
















Tables 
Structures For Actuator #1 Input Delay nk=1 
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After a set of candidate structures was obtained they were used to generate ARX models. 
The process of model selection and evaluation outlined in Figure 3.11 and described in 
sections 3 and 4 was applied to each candidate structure based ARX model. The most 
suitable candidate found was ARX [na=7 nb=6 nk=1] (see Figures 3.26 though 3.32). 
Table 3.3 lists the [na=7 nb=6 nk=1](or [7 6 1] ) structure as a manual selection. In the 
first iteration the [7 13 1] structure was selected. Upon repeating the selection process for 
actuator number two, the [7 6 1] structure was generated automatically by the algorithm 
(see appendix A) based on the choice of 13 parameters. The performance of the resulting 
model was judged superior and the [7 6 1] structure was applied to actuator number one 
where it also yielded the best ARX model (see Figure 3.32). However, two design trade- 
offs were made. First, the auto and cross correlation functions (see Figure 3.28) showed 
degraded performance, indicating there are still unaccounted input elements in the output. 
However, it does not affect the frequency response and other criteria used to validate the 
Structure. The second trade-off was the zero-pole plot (see Figure 3.31 ). The close 
proximity of two poles and two Zeros indicates that a lower order model would be 
desirable. However, this proved to not be the case. The use of a lower order model 


showed declining performance over the ARX model using the [7 6 1] structure. 
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Figure 3.26 
Frequency Response of Prefiltered ARX[7 6 1] (Actuator #1) 
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Plot of ARX MODEL [7 6 1] Output vs. Actual output “Data Prefiltered” 
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Figure 3.27 
Actuator #1 Output Vs. ARX[7 6 1] Model 
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Auto and Cross Correlation Functions of ARX[7 6 1} Model (Actuator #1) 
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Figure 3.29 
Residual Vs. Sample for ARX[7 6 1] (Actuator #1) 
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Figure 3.30 
Actual Output and ARX[7 6 1] Model Output Difference (Actuator #1) 


40 





Zero Pole Plot ARX Model [7 6 1] Model *Data Prefiltered* 
1.5 


-1 





Figure 3.31 
Zero-Pole Plot for ARX[7 6 1] Model (Actuator #1) 
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By increasing the delay the correlation functions improve but performance in all other 
respects is degraded. The 99% confidence limit may also be too constrictive. Confidence 
limits of 95% have been used [ref. 10 p.447] to constrain correlation functions. These 
factors governed the decision to accept the degraded correlation functions as a design 
trade off. In comparison with models based on various structures (Figure 3.32) this 


model is acceptable. Table 3.4 summarizes numerical results. 


nea 
Average Output Magnitude 0.1870 
Average Residual Magnitude 0.0064 


Average Output Difference Magnitude 0.0024 
Mean Square Fit 0.1562 


Table 3.4 
Prefiltered ARX[7 6 1] Numerical Results (Actuator #1) 





This process was repeated to find the structures for an ARX based model for actuator 
numbers two through six. This process involved the analysis of hundreds of 
validation/verification plots and numerical results. Only the validation plots of the 
accepted structure for each respective actuator were included in this thesis and are shown 


in appendix B. 


6. Model Structure Acceptance for Actuator Number Two 


The process outlined in section five was repeated for actuator number two. 
The plots resulting from the acceptance process are listed in appendix B (Figures A.1 
through A.7). Based on a delay of nk=1 from table 3.1 and resulting parameter number 
plot (see appendix B Figure A.1) the structures selected for evaluation are listed in table 


a. 
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Table 3.5 
Evaluated Model Structures for Actuator #2 










After evaluating the graphical and numerical results of the structures listed in table 3.5, 
the [7 6 1] structure was identified as the best structure. It provides excellent 
performance for a minimal number of zeros and poles. Figure 3.33 compares the [7 6 1] 


based model with the next best models. 
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Figure 3.33 
ARX Model Frequency Response Comparison 
Blue=Actuator, Red=[7 6 1], Magenta=[11 14 1], Black=[11 11 1] 
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Table 3.6 
Numerical Results for Actuator #2 


ie Model Structure Acceptance for Actuator Number Three 


The plots resulting from the acceptance process are listed in appendix B 
(Figures B.1 through B.7). Based on a delay of nk=1 from table 3.1 and resulting 
parameter number plot (see appendix B Figure B.1) the structures selected for evaluation 


are listed in table 3.7. 


Number of Parameters Structure [na nb nk] 


Table 3.7 
Evaluated Model Structures for Actuator #3 














A comparison of the selected structure, [7 6 1], and the next best models is shown in 
Figure 3.34. Of particular interest, the best structures were those resulting from non- 
filtered data. A possible explanation is the absence of high frequency disturbances which 


filtering is designed to remove. Also of note is the fact that during excitation of actuator 
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number three the noise emitted was distinctly different from the other actuators when 


excited. This could be the result of slight manufacturing differences. 
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ARX Model Frequency Response Comparison 
Blue=Actuator#3, Red=[7 6 1], Magenta=[15 15 1], Black=[7 13 1] 
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Average Output Magnitude 0.2298 


Average Residual Magnitude 0.0613 












Percentage of Average Output 26% 


Average Output Difference Magnitude 0.0037 
Mean Square Fit 0.2061 


Table 3.8 
Numerical Results For Actuator #3 





8. Model Structure Acceptance for Actuator Number Four 


The plots resulting from the acceptance process are listed in appendix B 
(Figures C.1 through C.7). Based on a delay of nk=1 from table 3.1 and resulting 
parameter number plot (see appendix B Figure C.1) the structures selected for evaluation 


are listed in table 3.9. 


Number of Parameters Structure [na nb nk] 






S12 2 MDL Based 


Se 
ie 


Table 3.9 
Evaluated Model Structures 










Similar to the case of actuator number one, the algorithm selected a delay of two. 


However, the best results were achieved by using a lower delay. The penalty paid was 
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degraded auto and cross correlation functions. This was again taken as a design trade off. 
The ARX structure [7a=7 nb=6 nk=1] again proved to be the best choice (see Figure 
3.35). Numerical results are contained in table 3.10. 
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Figure 3.35 
ARX Model Frequency Response Comparison 
Blue=Actuator#4, Red=[7 6 1], Magenta=[11 13 1], Black=[11 11 1] 
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Average Output Magnitude 0.3052 
Average Residual Magnitude 0.0066 


Average Output Difference Magnitude 0.0036 
Mean Square Fit 0.2662 


Table 3.10 
Numerical Results for Actuator #4 





a Model Structure Acceptance for Actuator Number Five 


The plots resulting from the acceptance process are listed in appendix B 
(Figures D.1 through D.7). Based on a delay of nk=1 from table 3.1] and resulting 
parameter number plot (see appendix B Figure D.1) the structures selected for evaluation 


are listed in table 3.12. 
















Table 3.11 
Evaluated Model Structures 


The validation of actuator number five is virtually the same as actuator number four (see 


appendix B Figures C.2 through C.7, D.2 through D.7 and tables 3.10 and 3.12). 
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Figure 3.36 
ARX Model Frequency Response Comparison 
Blue=Actuator #5, Red=[7 6 1], Magenta=[15 9 1], Black=[13 15 1] 
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Table 3.12 
Numerical Results for Actuator #5 








10. Model Structure Acceptance for Actuator Number Six 


The plots resulting from the acceptance process are listed in appendix B 
(Figures E.1 through E.7). Based on a delay of nk=1 from table 3.1 and resulting 
parameter number plot (see appendix B Figure E.1) the structures selected for evaluation 


are listed in table 3.13. 


Number of Parameters Structure [na nb nk] 
30 | 15151 ADC/MDL Based 


: 
7 





Table 3.13 
Evaluated Model Structures 


The comparison of the best model structures is shown in Figure 3.37 . The [7 6 1] 
structure selected for the other five actuators is the best choice for actuator number six. 


Numerical results are listed in table 3.14. 
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Figure! 3.37 
ARX Model Frequency Response Comparison 
Blue=Actuator #6, Red=[7 6 1], Magenta=[15 15 1], Black=[12 12 1] 


a3 





Average Output Magnitude 0.1551 
Average Residual Magnitude 0.0055 


Average Output Difference Magnitude 0.0023 
Mean Square Fit 0.1316 


Table 3.14 
Numerical Results for Actuator #6 
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D. APPLICATION OF SELECTED STRUCTURE 


The selection of the same basic model structure na=7 nb=6 nk=1 for all six 
actuators is an indication that the plant dynamics of the actuators can be assumed to be 
identical. With the exception of actuator number three, all received the best results from 
prefiltered data. The design trade-offs were the degraded auto and cross correlation 
functions and the relative proximity of two poles and two zeros. This basic structure was 
used to develop more accurate ARMAX and Box-Jenkins models. Although the Box- 
Jenkins is the most accurate parametric estimation method, it was anticipated that the 
ARMAX would give the best results based on the nature of disturbances and on when 
they enter into the system (see Figure 3.3). 

The application and evaluation process is the same as that used to evaluate 
selected structures on the ARX model. A second order noise model was used for both the 
ARMAX and Box-Jenkins models. The standard choices are 2" and 4" order models 
[Ref. 11]. Both were applied and the 2"° order yielded a better result. After the 
application of the ARMAX and Box-Jenkins based models the overall performance of the 
ARX and ARMAX models was decidedly better. The only exception was on actuator 
number two. This will be addressed in section two. Based on the results obtained thus 
far indicating the actuator plants are essentially the same and the excellent performance of 
the ARX and ARMAX models, only the evaluation plots for the ARMAX model are 
included and can be found in appendix C. Summary and comparison results are listed for 


each actuator. 


1. Actuator Number One 


A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.38. Validation plots for the ARMAX [7 6 2 1] model are included as 


Figures A.1 through A.6 in appendix C. Numerical results are contained in table 3.15. 
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Figure 3.38 
Model Frequency Response Comparison 
Blue=Actuator #1, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure ARX[7 6 1] ARMAX[7 6 2 1] BJ(6 2 27 1} 
_ (Data Prefiltered) (Data Prefiltered) (Data Prefiltered) 


Table 3.15 
Numerical Results for Actuator #1 





Based upon the above results the ARX[7 6 1] model was determined to best fit the 
performance of actuator number one. The choice between the ARX and ARMAX Is a 
close one. This indicates the assumption that the disturbance or noise entered the system 
early and is accurately modeled by the dynamics (poles) of the plant. The addition of the 
C(q,8) polynomial is not necessary. The improved performance of the ARX and 
ARMAX models over the Box-Jenkins further supports the assumption regarding the 


nature of the disturbance. 


2. Actuator Number Two 


A graphical comparison of ARX, ARMAX and Box-Jenkins models 1s 
shown in Figure 3.39. Validation plots for the ARMAX [7 6 2 1] model are included as 


Figures B.1 through B.6 in appendix C. Numerical results are contained in table 3.16. 
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Figure 3.39 
Model Frequency Response Comparison 
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Measure ARX[7 6 1] ARMAX[7 62 1] BJ[6 22 7 1) 
a. (Data Prefiltered) (Data Prefiltered) 


Table 3.16 
Numerical Results for Actuator #2 








The Box-Jenkins model showed excellent frequency response correlation, auto and cross 
correlation functions and some improved numerical results. However, due to the close 
performance of the ARMAX and Box-Jenkins models, the large uncertainty ellipses in 


the zero-pole plot (see appendix C) the ARMAX model was selected. 


3: Actuator Number Three 


A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.40. Validation plots for the ARMAX [7 6 2 1] model are included as 


Figures C.1 through C.6 in appendix C. Numerical results are contained in table 3.17. 
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Figure 3.40 


Model Frequency Response Comparison 
Blue=Actuator #3, Reda=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure ARX[7 6 1] ARMAX[(7 6 2 1] BJ[6 227 1] 
— (Data not filtered) (Data Prefiltered) 


Table 3.17 
Numerical Results for Actuator #3 








Although the non-filtered ARX model was chosen over the prefiltered ARX model, it 
appears that the performance of the prefiltered ARMAX is much better. This can be 
explained by possible manufacturing variation in actuator number three which requires 


the use of the C(q,8) polynomial to model the disturbance. 


4. Actuator Number Four 


A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.41. Validation plots for the ARMAX [7 6 2 1] model are included as 


Figures D.1 through D.6 in appendix C. Numerical results are contained in table 3.18. 
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Figure 3.41 
Model Frequency Response Comparison 
Blue=Actuator #4, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 
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Measure ARX[7 6 1] ARMAX[7 62 1] BJ[6 227 1] 
a (Data Prefiltered) (Data Not Prefiltered) 
Ave. ee rt Magaitnde Magnitude p 03052 p 03052 os p 03052 


Table 3.18 
Numerical Results for Actuator #4 





This case appears similar to actuator number one where the disturbance 1s accurately 


characterized by the system pole model without the need for the C(q,8) polynomial. 


S Actuator Number Five 


A graphical comparison of ARX, ARMAX and Box-Jenkins models is 
shown in Figure 3.42. Validation plots for the ARMAX [7 6 2 1] model are included as 


Figures E.1 through E.6 in appendix C. Numerical results are contained in table 3.19. 
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Figure 3.42 
Model Frequency Response Comparison 
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Measure ARX[7 6 1] ARMAX[7 62 1] BJ[6227 1] 
(Data filtered) (Data Prefiltered) (Data Pre filtered) 
Ave. | Ave. Output Magnitude | Magnitude p 03070, 0.3070 0.3070 


a 0.0035 0.0151 


Ave. Output Difference 0.0032 0.0033 0.0034 
Mean Square Fit 0.2830 0.2688 0.2698 


Table 3.19 
Numerical Results for Actuator #5 








Ave. — Mag. 








6. Actuator Number Six 


A graphical comparison of ARX, ARMAX and Box-Jenkins models 1s 
shown in Figure 3.43. Validation plots for the ARMAX [7 6 2 1] model are included as 


Figures F.1 through F.6 in appendix C. Numerical results are contained in table 3.20. 
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Figure 3.43 
Model Frequency Response Comparison 
Blue=Actuator #6, Red=ARMAX[7 6 2 1], Magenta=ARX[7 6 1], Black=BJ[6 2 2 7 1] 


66 













Measure ARX[7 6 1] ARMAX[7 62 1] BJ(6 227 1] 


Table 3.20 
Numerical Results for Actuator #6 








IE: SUMMARY OF RESULTS 


Table 3.21 presents a summary of the model selection process. 


er SE a a a ae a ee ee ee 
ARMAX 
[762 1] 


Model ARX ARMAX | ARMAX ARX ARMAX 
[Structure] [76 1] Wage li | (762 1) (76 1] [7621] 
0.1870 0.1906 0.2298 0.3052 0.3070 0.1551 
0.0064 0.003 0.0040 0.0066 0.0035 0.0035 
Mag. 
0.0024 0.0027 0.0040 0.0036 0.0033 0.0020 
Dive. 
0.1562 0.1705 0.2031 0.2662 0.2688 TEI 


Table 3.21 
Summary of Numerical Results 








F. POLES AND ZEROS 


The ultimate goal in identifying the plant of the platform is to obtain the poles and 
zeros of the plant transfer function which can be used in design of an active vibration 
control law for the platform. Tables 3.22 and 3.23 present the poles and zeros extracted 


from the parameter estimation algorithm listed in appendix A. 
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0.1252+0.91991 0.1080+0.96671 0.1198+0.92041 
0.1252-0.91991 0.1080-0.96671 0.1198-0.92041 


1.0137 0722 1.0204 
0.5371+0.61671 0.5887+0.64861 0.4925+0.62061 


0.5371-0.61671 0.5887-0.64861 0.4925-0.62061 


(a) Actuators #1 through #3 


eC 


(b) Actuators #4 through #6 


Table 3.22 
Actuator Transfer Function Poles 
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a 


(a) Actuators #1through #3 
I ee ee 
0.6905 + 0.66661 0.1657 + 0.90001 0.1138 + 0.8029) 





























(b) Actuators #4 through #6 


Table 3.23 
Actuator Transfer Function Zeros 


A review of the poles and zeros shows that zeros are very similar, especially among 
actuators one, two, three and six. The magnitude of the real parts of zeros for actuators 
four and five appear to vary more than the others. This is also true for the poles. This 
was unexpected and no explanation is offered. Any expected difference would have been 
attributed to actuator number three based on previous results and analysis. However, 


actuator number three is in relative concurrence with the other actuators. 
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G. TRANSFER FUNCTIONS 


The actuator transfer functions are as follows; 


ACTUATOR Number One 


13.5362S° —31.6493S° + 42.5317S* —39.4886S° + 22.8012S* —7.9085S 
S? —3.8132S° + 7.7793S8° —10.35325* + 9.7102S° —6.3907S7 + 2.73365 — 0.6064 


ACTUATOR Number Two 


15.1709S° —36.6464S° + 51.4596S* —49.9296S? + 30.8599S? —11.25845 
S’? —3.7010S° + 7.5692S° —10.2511S* +9.9482S° — 6.87628? + 3.13325 —0.7594 


ACTUATOR Number Three 


21.6875S° —48.6904S° + 64.51985* —59.84855S° + 33.8370S7 —11.9673S 
S? —3.2342S° +6.0679S° — 7.5022S* +6.6380S° — 4.11758” + 1.64545 — 0.3361 


ACTUATOR Number Four 


13.4875S° —30.7993S° + 40.4032,5* —36.9017S° + 21.2236S? — 7.5489S 
S7 —3.9270S° +8.0953S° —10.8076S* +10.1338,5° — 6.634887 + 2.80155 — 0.6022 


ACTUATOR Number Five 


8.51315° + 3.9823S8° — 9.42375* — 4.22425? + 0.930025" —0.1950S 
S’ —2.0359S° +1.7438° —0.4824S° +0.0175S° —0.27736- +0257 750 so 


ACTUATOR Number Six 


9.6768S° —11.5692S° + 11.0591S8* —9.1858S° + 2.52278? —2.5671S 
S' —2.7542S8° +4.5011S° —4.7744S* + 3.6582S° —1.9692S" +0.68115 —0.1253 


There is greater concurrence among actuators one, two, three and six. However, overall 
comparison of zeros, poles and transfer functions indicates that all six plants are very 


similar. 
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IV. ACTUATOR COUPLING AND VIBRATION CONTROL 
APPLICATION 


A. ACTUATOR COUPLING 


With the dynamic models constructed in Chapter Three it is now possible to model 
and investigate the interaction between actuators to see if indeed they can be controlled 
independently. The basic ARX method has been proved to be quite accurate in modeling 
the plant transfer function of each actuator. To determine the presence of actuator 
coupling each actuator was excited with 50mV pink noise and the output response of all 
six actuators measured. The input-output data gathered was introduced to the ARX[7 6 
|] model and the resulting model frequency response was plotted. Figures 4 | through 4.5 


show the results from the excitation of actuator number one. 
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Figure 4. | 
Frequency Response Actuator #1 Vs. Actuator #2 
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Figure 4.2 
Frequency Response Actuator #1 Vs. Actuator #3 
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Figure 4.3 
Frequency Response Actuator #1 Vs. Actuator #4 
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Frequency Response Actuator #1 Vs. Actuator #5 
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In Figures 4.1 through 4.4 there is discernible correlation between the frequency 
response of actuator number one and the frequency response of actuators two through five 
to the excitation of actuator number one. Of particular interest is the response at the 
resonance frequency of 1.4 kHz. On average there is approximately a 10 dB difference at 
1.4 kHz for these actuators. In the frequencies below 1kHz the response is in the 
negative dB range. This would indicate that the passive stage is performing well in 
eliminating the low frequency disturbances. Although coupling is present, it is relatively 
low in comparison to the response of actuator one. However this is not the case with 
actuator number six. As in actuators two through five the response below 1 k is small. 
Above IkHz, particularly at 1.4 kHz the response is close to a mirror image of the 
response of actuator number one. This suggests that there is significant high frequency 


coupling 


B. CONTROL APPLICATION AND COUPLING VERIFICATION 


Given the results of the previous section it is necessary to verify on the UQP. 
This was accomplished by applying a control law to the platform and checking for any 
evidence of coupling. The first evaluation conducted consisted of controlling an 
individual actuator and the second applied control to all six actuators simultaneously. 
There are numerous control candidates for use in AVC. For random broad band 
noise/vibrations, feedback control is preferable. Feed forward controllers are good 
candidates when a specific frequency of the disturbance is known. Local (decoupled) 
force feedback [Ref. 5], adaptive vibration control using two least mean square filters 
[Ref. 4], and absolute velocity feedback [Ref. 6] are examples that have been used in 
AVC. The latter of the three also used a geophone for sensing. The control law currently 
programmed for the system is a lead-lag compensator. This type of compensator is a 
variant of derivative and integral control. Lead compensation has the effect of lowering 
rise time and decreasing overshoot while lag compensation is used mainly to gain steady 
State accuracy of the system[Ref. 7]. A lead-lag controller was provided by CSA 


Engineering and was used to ensure the operation of each individual strut. This control 
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law was modified to apply the same lead-lag compensator to all six actuators 


simultaneously based on the initial assumption of independently controllable actuators. 


The results previously obtained indicate that this cannot be done. 


Compensator Transfer Function 


1. 


A Bode Plot of the magnitude and phase of the compensator (See Figure 


4.6) shows an infinite gain margin and a phase margin of approximately 60 degrees. 
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Figure 4.6 
Compensator Bode Plot 
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A bode plot of the open loop tranfer function is provided in Figure 4.7. 
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Figure 4.7 
Open Loop Transfer Function (Compensator* Actuator) 


Finally, a magnitude plot in terms of dB of closed to open loop is shown in Figure 4.8. 


From this plot it is observed that the best performance for this compensator will be 


between 12 Hz to 48 Hz with maximum reduction of 32dB at 25 Hz. A reduction of 


This controller was designed with an 


approximately 25 dB can be seen at 42 Hz. 


approximate bandwidth of 100 Hz. The spike at approximately 10 Hz is a faulty value 


resulting from system noise or an error in data collection. 
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Figure 4.8 


Closed Loop/Open Loop 


Control Application 


Pa 


applied to actuator 


The controller described in section one was first 


number two and then the .modified controller was applied to all six actuators 


simultaneously. Table 4.1 details the equipment preparation and settings. 


Ue 


Equipment Setting 


HP 33120A Function Generator Waveform Sinusoid 
Frequency 42 Hz 
Amplitude 500mv peak- 
to-peak 
KEPCO Bipolar Operational Power Current Control Mode 
Supply/Amplifier 





Table 4.1 
Equipment Settings 


3, Experiment 


The first task is to look at the platform in open loop mode. From this the 
local environment can be assessed. The plots in Figures 4.9 and 4.10 show what is 
sensed by the geophone on actuator two with the 42 Hz disturbance and no control (open 


loop). 
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Figure 4.9 
Ambient PSD (0-5kHz) 
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Open Loop-42Hz Disturbance 
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Figure 4.10 
Ambient PSD (0-200Hz) 
Of note in Figure 4.9 is the 1.4 kHz resonance frequency. This is the result of power 
being applied to the actuator from the HVPS. At approximately 30 Hz there is an 
environmental disturbance of unknown origin. In the open loop mode the 42Hz 
disturbance is clearly visible (see Figure 4.10). Figure 4.11 shows the power spectral 
density (PSD) of the sensor output for actuator number two for open and closed loop 
operations (both single and all six actuators). With only power applied to the 
piezoceramic actuator (no control), excitation of the resonance frequency is visible as the 
blue trace in Figure 4.11(a). Closing the control loop for single and multi-actuator 
control increases the excitation. Figure 4.12 shows the sensor output for actuators two, 
three, and six with control applied to actuator number two. Actuator number three shares 
the same base node with actuator number two. The coupling between actuators two and 
three can be seen in Figure 4.12(a). The application of control excites the resonance of 
all three actuators, but the interaction between actuators two and three is the most 


pronounced. 
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(b) PSD (0-200Hz) Sensed By Geophone on Actuator #2 


Figure 4.11 
Open Vs. Closed Loop Performance PSD 
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(b) PSD (0-200Hz) Sensor Outputs 


Figure 4.12 
Actuator #2 Closed Loop Performance PSD 
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Vv. SUMMARY 


A. CONCLUSIONS 


The process of system identification provided valuable insight into the UQP. By 
treating the UQP as a “Black Box” and using ready-made models the complex system 
dynamics were overcome and very accurate models of the system were developed. Some 
possible reasons for differences between the models and the actual system can be 
explained by non-removal of faulty values. Due to the number of samples used to build 
and validate the models, hand smoothing of the data was not practical. Based on the 
average residuals for each actuator being less than 4% of the average model output, 
outliers did not pose a significant problem. In comparing the validation data (numerical 
and graphical), actual frequency responses, and zeros, poles and transfer functions there 
are distinct similarities. The fact that all six are described by seven poles and six zeros 
was anticipated. However, some variation among actuators was expected. There is 
strong correlation of actuators one, two, three and six are. Although the numerical values 
for the zeros, poles and transfer function of actuators four and five did not correlate with 
the other actuators as strongly, they still retained the same pole-zero structure. Based on 
this it can be stated that for the purposes of characterizing the system and control design, 
all six actuators are identical. 

The sensitivity of the geophones resulting in the introduction of external 
disturbances could be a source of inaccuracy, however, they are considered small in 
comparison to the disturbances generated by the other actuators. The emergence of ARX 
and ARMAX models as the most accurate models validates the assumption that these 
disturbances to the system entered the process early on. The disturbance created by the 
excitation of an actuator in turn causes excitation of other actuators, especially for 
actuators that share the same aluminum base in the region above 200 Hz. Since the 
disturbance is caused by another actuator, the choice of modeling the noise with the 
system dynamics (poles) was an accurate one. The interaction or coupling among 
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actuators, particularly between base adjoining actuators is distinct. Although the passive 
damping mechanism works well in eliminating low frequency interaction, high frequency 
disturbances are readily transmitted between actuators. For actuators not sharing a 
common base node this interaction is minimal and could possibly be ignored. However, 
for base node sharing elements this problem is significant. 

The lead-lag compensator worked well for the single actuator in the 1Hz to 100Hz 
range, providing over 20dB damping of the 42Hz disturbance. However, above 200 Hz 
the controller exacerbates the disturbances. This problem is amplified in the case where 
all six actuators are controlled simultaneously. The interaction between actuators 
compounds the high frequency degraded performance of the controller. This causes it to 
go unstable. The high frequency interaction between base sharing nodes must be 


addressed before attempting to control all six actuators simultaneously. 


B. RECOMMENDATIONS FOR FUTURE WORK 


With accurate plant models of the actuators the next logical step is control law 
design. With knowledge of the system it is possible to design controllers using either 
classical or modern techniques. The main concern in the design is the high frequency 
interaction in the range of 1kHz to 1.5 kHz. The revelation of actuator interaction does 
not necessarily mean that it is impossible to control each actuator independently. The 
dynamic compensator applied to an individual actuator worked well below 100Hz but did 
not roll off sufficiently. This resulted in high frequency instability in control of all six 
actuators. With knowledge of the plant transfer function it is possible to design a 
dynamic compensator capable of controlling each actuator independently. However, the 
design must be such that it does not excite the 1.4 kHz resonance. On the other hand, it 
may be necessary to treat the UQP as a multiple input/multiple output (MIMO) system 
and a modern control design is therefore required. Another possible remedy to the 
interaction between base node sharing elements is to change the material of the 
attachment node. The use of aluminum acts as an excellent conductor of vibration from 
one actuator to another. By using vibration absorbing material such as _ high stiffness 


rubber it may be possible to reduce the interaction sufficiently such that each actuator 
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could be controlled independently using single input/single output (SISO) classical design 
techniques. Replacing the upper node material could also reduce the overall transmission 


of unwanted vibration, not only from interaction but external sources as well. 
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APPENDIX A. COMPUTER PROGRAMS 


All computer codes included in this thesis were written by Naval Postgraduate 
School Space Research and Design Center unless otherwise noted. The code included in 
this appendix 1s organized along lines of the thesis. 


A. SYSTEM IDENTIFICATION AND MODEL CONSTRUCTION 

% sysident.m : MAIN PROGRAM 

% Written by: George D. Beavers June 1997 

% This MATLAB “.m” file reads in the user requested actuator input output data 
% from the “4h*.mat”’ files. “*’’ represents the actuator that is being excited by the 
% SOmvV “pink noise’. It calls the user written functions 

% “arxmod.m”, “armaxmod.m” and “bjmod.m” to calculate models based on the 
% ARX, ARMAX and Box-Jenkins parameter estimation methods respectively 

% ‘The zeros, poles and transfer functions are then extracted from the models 

% returned by the user written functions. 

% Each “4h*.mat” file has 7 rows of data 

%o rows 1-6 = geophone output from the respective actuator (row 1=actuator #1) 
% row 7 is the 50mv “pink noise” input 

% Information on the MATLAB SITB functions can be found in reference[12] 


clear 

dum='y'; 

while (dum=='y’), 

j=input(‘Enter strut number >> '); 
disp([‘load 4h’ int2str(j)]}) 


eval([‘load ' '4h' int2str(j) ]) % Load input output data % 
z2=[trace_y(j,1:10000)' trace_y(7,1:10000)']; %oPull out requested channel % 
z2=dtrend(z2); %Remove trends & adjust levels% 


z3=[trace_y(j, 10001:20000)' trace_y(7, 10001 :20000)']; 
z3=dtrend(z3); 


disp(['Choose Structure Selection Method’]) 


disp({'l1==Automatic’]) % Use algorithm to determine % 
disp({'2==Manual')) %[na,nb, nk] orinputthem % 
temp=input('>>> '); % manually % 


disp({'Choose a Parametric Estimation Model']) 
Ainpaiiee——o5 “armax =2'; 9 bj =3')) 
che=input('>>> ') 


% The following subroutine is used to determine input delay and number of parameters % 
% and structure na,nb based on that delay for an ARX based model% 
if temp==1 

V=arxstruct(z2,z3,struc(2,2,1:5)); 
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[nn, VmJ=selstruc(V, AIC’) 
inp=input(‘enter delay’) 
V=arxstruct(z2,z3,struc(1:15,1:15,inp)); 
nnAsselstruct(V,'AIC’) 
nnMe=selstruct(V, MDL’) 
Figure(1) 
nns=selstruct(V) % Generate parameter Vs. Loss fen% 
% plot Jo 
% Determine which model to construct then call the function% 
% Using the automatically generated ARX based Structure % 
ir ch¢==— 
arxmod(j,nns,z2,z3); 
elseifiche——2 
nc=input(‘Enter Order of Noise model nc >>"); 
itr=input(‘Enter Maximum Number of Iterations >>’); 
narm=[nns(1,1) nns(1,2) nc nns(1,3)]; 
armaxmod(j,narm,itr,z2,z3) 


else 
on=input('Enter Order of Noise model [nc nd] >>") 
itr=input(‘Enter Maximum Number of Iterations >>") 
nbj=[nns(1,2) on nns(1,1) nns(1,3)] 
bjmodQG,nbj,itr,z2,z3) 

end 


% Generate models based on manually entered structure % 
else 


i chce——) 
nns=input(‘Enter Structure [na nb nk] >>’); 
[arxm,arxmft]=arxmod]1(j,nns,z2,z3); 
[zepol ,k1]=th2zp(arxm) % Extracts zeros-poles 
[num 1,den] ]=th2tf(arxm) % Extracts transfer function 
[zepolf,k1f]=th2zp(arxmft) % Extracts zeros-poles (filtered data) 
[num1f,den1f]=th2tf(arxmft) 9% Extracts tranfer function(filtered data) 
present(arxm) % Extracts parameters 
present(arxmft) % Extracts parameters(filtered data) 
elseil che——2 
narm=input(‘Enter Structure [na nb nc nk] >>"); 
i1tr=input(‘Enter Maximum Number of Iterations >>') 
[armxm,armxmft]=armaxmdf(j,narm, itr,z2,z3); 
[zepo2,k2]=th2zp(armxm) % Extracts zeros-poles 
[num2,den2]=th2tf(armxm) 9% Extracts tranfer function 
[zepo2f,k2f]=th2zp(armxmft) % Extracts zeros-poles (filtered data) 
[num2f,den2f]=th2tf(armxmft) % Extracts tranfer function(filtered data) 
present(armxm) % Extracts parameters 
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present(armxmift) % Extracts parameters(filtered data) 


else 
nbj=input(‘Enter Structure [nb nc nd nf nk] >>") 
itr=input(‘Enter Maximum Number of Iterations >>") 
[bjm,bjmft]=bjmdf(j,nbj,itr,z2,z3); 

end 


end 
dum=input(‘Continue ?>> ','s’); 
end 


To To To To To Yo To Yo To To Yo Yo Yo Yo To Yo To Yo Yo To Yo Yo To Vo Yo Yo Yo Lo To To To To To To To To To To Yo To Yo To Vo 


% arxmod.m % 
% ARX Model Construction and Analysis % 
%o George D. Beavers % 
Jo June 1997 %o 


To To To To To To To Yo To Yo Yo To To Fo To To Yo Fo To To To Yo Yo Yo To To Yo Yo To Vo Yo To To To To Yo Yo To Yo To Yo Yo Yo 
% This program generates a model based on the input-output data of asystem using % 
% the ARX parameter estimation method. It also generated a series of graphicaland % 
% numerical results in order to evaluate the suitability of the model Jo 

Jo To To To To To Yo To To To Yo To Yo Yo To To To To To To To Yo Yo Yo To To To Yo Yo To Yo To To Yo Yo To Yo Yo Yo Yo Yo Yo Yo 


function [arxm,arxmft]=arxmod(j,nns,z2,z3) 
Jo To To To To To Yo Yo Lo Lo To Qo Yo Fo Fo Yo Lo Do Yo Yo Yo Yo Yo Vo Yo Yo 


% arxm=ARX based model %o 
% arxmft=ARX based model (filtered data) % 
% j=actuator % 
% mnns=structure((na nb nc}) Jo 
% z2=input-output data for model construction % 
% Z3=input-ouput data for model evaluation Jo 


Zo To To To To To To Yo To To To To Yo Yo To Yo Yo Yo Yo To To Yo Yo Yo Yo Yo 


na=nns(1,1); 
NA=int2str(na); 
nb=nns(1,2); 
=int2str(nb); 
nk=nns(1,3); 
NKz=int2str(nk); 


arxm=arx(z2,nns); 
arxm=sett(arxm,0. le-3); % Set sample interval % 
garxm=th2ff(arxm); % Convert to frequency function% 


% This section loads separate data to compute the actual actuator frequency response% 


% and then plot it against the model’s frequency response Jo 
Figure(2) 

eval([‘load s',int2str(j),'m']) % 

X 1=0211x; 


¥ 1=20.*log10(abs(o0211)); 


semilogx(X1,Y1,'b’), 

hold on 

semilogx(garxm(:, 1)/(2*pi),20*log 10(garxm(:,2)),'r’); 

title({"Blue=Actuator#' num2str(j),,Red=ARX MODEL]',NA,' ',NB,' '.NK,']']) 
xlabel(Frequency (Hz)') 

ylabel(‘Magnitude(dB)') 

grid 

hold off 


Figure(3) 


eval([‘load s',int2str(j),'p']) 
X2=0211X; 
Y2=180*unwrap(angle(o211))/p1; 


semilogx(X2,Y2,'b') 

hold on 

semilogx(garxm(:,1)/(2*pi),(garxm(:,4)),'r'); 

title({‘Blue=Actuator#' num2str(j),', Red=ARX MODELT[',NA,' ',NB,'',NK,']']) 
xlabel('Frequency (Hz)’) 

ylabel('Phase (Deg)’) 

grid 

hold off 


% The above process is repeated using prefiltered data. A 10" order Butterworth filter % 
% is used. The passband is given in terms of fractions of the Nyquist frequency (0.52)% 
% This is a low pass filter with an upper bound of approximately 1.5kHz 

Jo 

zf=idfilt(z2,10,.52); %Prefilter model construction data% 


arxmft=arx(zf,nns); 
arxmft=sett(arxmift,O. le-3); 


garxmft=th2ff(arxmft); 


eval([‘load s',int2str(j),'m']) 
X1=0211x; 


O2 


Y 1=20.*log10(abs(0211)); 
Figure(4) 


semilogx(X1,Y1,'b’) 

hold on 

semilogx(garxmft(:,1)/(2*pi),20*log10(garxmft(:,2)),'r'); 
title([/Blue=Actuator#',num2str(j),', Red=ARX MODEL[',NA,' '|NB,' '|NK,'] *Data 
Prefiltered*']) 

xlabel(‘Frequency (Hz)’) 

ylabel(Magnitude(dB)’) 

grid 

hold off 


Figure(5) 


eval([‘load s',int2str(j),‘p']) 
~2=0211x; 
Y2=180*unwrap(angle(o2i1))/pi; 
semilogx(X2, Y2,'b') 


hold on 

semilogx(garxmft(:,1)/(2*pi),(garxmft(:,4)),'r'); 

title(['Blue=Actuator#' num2str(j),', Red=ARX MODEL[',NA,' ',NB,’ ',NK,'] *Data 
fretiitered* ’]) 

xlabel(‘Frequency (Hz)') 

ylabel(‘Phase (Deg)’) 

grid 

hold off 

% This section generates plots and numerical results used in evaluating the model % 
Figure(6) 

Yo Apply the input row of the validation data to the model and then compare with % 
ZY the actual output of the system corresponding to the validation data and plot results% 
[yh,fit]=compare(z3,arxm); 

axis([5000,5 100,-1.0,1.0}) 

title({'Plot of ARX MODEL [',NA,' ',NB,'',NK,'] Output vs. Actual output’]) 
text(5020,0.75,'Red=Model Output, Blue= Measured Output’) 

xlabel(‘Sample’) 

ylabel(‘Output Magnitude’) 

Figure(7) 


[yhf, fitf]=compare(zfv,arxmft); %Prefiltered data% 


title({'Plot of ARX MODEL [',NA,’ ',NB,' ',NK,'] Output vs. Actual output *Data 
Prefiltered*']) 
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axis({5000,5100,-1.0,1.0]) 

text(5020,0.75,'Red=Model Output, Blue= Measured Output’) 
xlabel(’Sample’) 

ylabel(‘Output Magnitude’) 


% Calculate the residuals and plot the auto and cross correlation functions % 


Figure(8) 

el=resid(z3,arxm); % Non-filtered data 
Figure(9) 

e2=resid(zfv,arxmft); % Prefiltered data 


% Plot the residuals % 

Figure(10) 

plot(el,'m’),grid 

title({"Residual of ARX Model [',NA,' ',NB,' '",NK,']']) 
ylabel('Residual’) 

xlabel('Sample’) 

axis({5000,5 100,-0.1,0.1]) 

Figure(11) 

plot(e2,'m’),grid 

title({'Residual of ARX Model [',NA,' ',NB,' ',NK,'] *Data Prefiltered*']) 
axis([5000,5 100,-0.05,0.05]) 

xlabel(‘Sample’) 

ylabel(‘Residua!’) 


% Calculate Numerical Results% 
el pos=abs(e1); 

e2pos=abs(e2); 
ypos=abs(z3(:,1)); 


ave=mean(e 1 pos) % Average residual 
aveft=mean(e2pos) % Average residual 
outave=mean(ypos) % Average actual output 

fit % Mean Square Fit calculated above 
fitf %Mean Square Fit calculated above 


% Calculate and plot the actual output-model output 

Z30ut=z3(:, 1); 

for 1=5000:5 100; 
difft(ij=z3 out(i)-yhf(i); 

end 

Z50Ut=73(6, 1); 

for j=5000:5 100; 
diff(j)=z3 out(qy)-yhq); 

end 

Figure(14) 

plot(diff,'r’) 
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axis({5000,5 100,-0.5,0.5]) 

title({'Actua] Output - ARX Model [',NA,' ',NB,' '.NK,'] Model']) 

xlabel(‘Sample’) 

Figure(15) 

plot(difft,'r’) 

ax1s({5000,5100,-0.5,0.5]) 

title({'Actual Output - ARX Model [',NA,’ ',NB,' '|NK,'] Model *Data Prefiltered*'}) 
xlabel(‘Sample’) 

avediff=mean(abs(diff)) YoAve. difference between actual and model outputs 
avedft=mean(abs(difft)) Yo Ave. difference between actual and model outputs (filter) 


% Generate Zero-Pole plots 

Figure(12) 

zpplot1(th2zp(arxm),3,{],1.2) 

title({"Zero Pole Plot ARX Model [',NA,’ ',NB,' 'J,NK,'] Model']) 

Figure(13) 

zpplotl(th2zp(arxmft),3,[],1.2) 

title({'Zero Pole Plot ARX Model [',NA,' ',NB,’ ',NK,'] Model *Data Prefiltered*']) 


JoTo Fo To To Po Po To To To To Fo Fo Po To To To To To To To To Fo Fo To To To To To Lo Fo To To Fo Mo To Yo Vo Fo Vo ToVo Fo 


% armaxmod.m Jo 
% ARMAX Model Construction and Analysis Jo 
Yo George D. Beavers Jo 
% June 1997 J 


Jo To To To To To Yo To To Yo To Yo Yo Yo Yo To To To To To To Po To To Yo Yo To Yo To To To To Yo Yo Yo Yo To To Fo Po Po Po Fo 
% This program generates a model based on the input-output data of asystemusing % 
% the ARMAX parameter estimation method. It also generated a series of graphical % 
% and numerical results in order to evaluate the suitability of the model J 
Jo To To To Yo To To No To To Yo To Fo To Yo To To To Fo Po To To Fo Fo Yo Po To Po Yo Po To Po To Po Fo Po Yo Po Yo Yo Po Yo Yo 


function [arxm,arxmft]=armaxmod(j,nns,itr,z2,Z3) 
Jo To To To To To To Po To Po Po To To Fo Yo To To To To Po To Fo Fo Fo Po%o 


Yo arxm=ARMAX based model Jo 
% arxmft=ARMAX based model (filtered data) Jo 
% j=actuator Jo 
% mnns=structure({na nb nc}]) Jo 
%o z2=input-output data for model construction Zo 
% Z3=input-ouput data for model evaluation Jo 


To To To To To To Yo Fo To To To To Yo To To To To To Yo To Yo Yo Fo To To Vo 


na=nns(1,1); 
NA=int2str(na); 
nb=nns(1,2); 
NB=int2str(nb); 
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nc=nns(1,3); 
NC=int2str(nc); 
nk=nns(1,4); 
NKe=int2str(nk); 


% This section loads separate data to compute the actual actuator frequency response% 
% and then plot it against the model’s frequency response Jo 


arxm=armax(Z2,nns, itr); 
arxm=sett(arxm,0. 1le-3); % Set sampling interval Jo 
garxm=th2ff(arxm); % Convert to frequency function% 


eval([‘load s',int2str(j),'m']) % Load actual actuator data (magnitude)% 
X 1=0211x; 

Y 1=20.*log10(abs(0211)); 

Figure(2) 

semilogx(X1,Y1,'b)), 

hold on 

semilogx(garxm(:,1)/(2*pi),20*log10(garxm(:,2)),'r'); 

title(['Blue=Actuator#' ,num2str(j), Red=ARMAX MODEL[',NA,'',NB,' ',NC,' ',NK,']']) 
xlabel(‘Frequency (Hz)’) 

ylabel(Magnitude(dB)’) 

grid 

hold off 


eval([‘load s',int2strG),'p']) % Load actual actuator data (phase)% 
X2=0211x; 

Y2=180*unwrap(angle(o211))/pi; 

Figure(3) 

semilogx(X2, Y2,'b’) 

hold on 

semilogx(garxm(:,1)/(2*pi),(garxm(:,4)), 1’); 
title({'Blue=Actuator#' ,num?2str(j), Red=ARMAX MODEL[',NA,'',NB,' ',NC,' ',NK,'J']) 
xlabel(‘Frequency (Hz)') 

ylabel(‘Phase (Deg)’) 

grid 

hold off 


Op #8 KEK KEKE ARMAK MODEL (DATA PREFILTERED ***% *## #4 HH EE EEK KKK Of, 


% The above process is repeated using prefiltered data. A 10" order Butterworth filter % 
% is used. The passband is given in terms of fractions of the Nyquist frequency (0.52)% 
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7 This is a low pass filter with an upper bound of approximately 1.5kHz 
%o 

7i=idfilt(z2,10,.52); % Prefilter model construction data % 
arxmft=armax(zf,nns, itr); 

arxmft=sett(arxmft,0. le-3); 

garxmft=th2ff(arxmft); 


eval([‘load s',int2str(j),'m']) 
peli=o211x; 
Y 1=20.*log10(abs(02i1)); 


Figure(4) 

semilogx(X1,Y1,'m’) 

hold on 

semilogx(garxmft(:, 1)/(2*pi),20*log 10(garxmft(:,2)),'2'); 

title({(‘Blue=Actuator#' num2str(j),,Red=ARMAX MODEL[',NA,' ',|NB,' NC,’ '"|NK,'] 
*Data Prefiltered*']) 

xlabel('Frequency (Hz)') 

ylabel(‘Magnitude(dB)’) 

grid 

hold off 


eval([‘load s',int2str(j),'p']) 
pe7=0211x: 
Y2=180*unwrap(angle(o0211))/pi; 


Figure(5) 

semilogx(X2, Y2,'b') 

hold on 

semilogx(garxmft(:,1)/(2*p1),(garxmft(:,4)),'r'); 
title({'Blue=Actuator#' num ?2str(j),,Red=ARMAX MODEL[',NA,'',NB,'',NC,' 'JNK,'] 
*Data Prefiltered*'}) 

xlabel(‘Frequency (Hz)') 


ylabel(‘Phase (Deg)’) 

grid 

hold off 

To ¥ KKK KA KAKKEKEEEEKGR APHICAL AND NUMERICAL ANALY SIS 42 F2F ts tttt@ 
% This section generates plots and numerical results used in evaluating the model % 
% Apply the input row of the validation data to the model and then compare with Jo 


% the actual output of the system corresponding to the validation data and plot results 
%o 

Figure(6) 

[yh,fit]=compare 1(z3,arxm); 

axis({5000,5100,-1.0,1.0}) 


7 


title({‘Plot of ARMAX MODEL[',NA,’ ',NB,'',NC,' ',NK,'] Output vs. Actual output')) 
text(5020,0.75,'Red:Model Output, Green: Measured Output’) 

grid 

Figure(7) 

[yhf,fitf]=compare | (zfv,arxmft); 

title([‘Plot of ARMAX MODEL|',NA,’',NB,’',NC,' ',NK,"] Output vs. Actual output 
*Data Prefiltered*'}) 

axis({5000,5100,-1.0,1.0]) 

xlabel(‘Sample’) 

grid 


% Calculate residuals and generate and plot auto & cross correlation functions % 
Figure(8) 

el=resid1(z3,arxm); 

Figure(9) 

e2=resid1(zfv,arxmft); 


% Plot Residuals % 

Figure(10) 

plot(el,'m’),grid 

title(["Residual of ARMAX MODEL|’,NA,' ',NB,' ',NC,' ',NK,']']) 
ylabel('Residual’) 

ax1s({[5000,5 100,-0.1,0.1]) 

Figure(1 1) 

plot(e2,'m’'), grid 

title({"Residual of ARMAX MODELT|',NA,’ ',NB,'',NC,’ 'JNK,'] *Data Prefiltered*']) 
ax1s({5000,5 100,-0.05,0.05]) 

xlabel(‘Sample’) 

ylabel('Residual’) 


% Calculate Numerical Results% 
el pos=abs(e1); 

e2pos=abs(e2); 
ypos=abs(z3(:,1)); 


ave=mean(e 1 pos) % Average residual 
aveft=mean(e2pos) % Average residual (Prefiltered) 
outave=mean(ypos) % Average actual output 

fit % Mean Square Fit calculated above 
fitf % Mean Square Fit (Prefiltered) 


% Calculate actual output minus model output an plot results % 


Z30ut=Z3(:,1); 
for 1=5000:5 100; 
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difft(ij=z3 out(1)-yhf(i); 
end 
Zeout=73(:,1); 
for j=5000:5 100; 
diff(j)=z3 out(j)-yhQ); 
end 
Figure(14) 
plot(diff,'r’) 
axis({5000,5 100,-0.5,0.5]) 
title([{'Actual Output - ARMAX[',NA,'',NB,' ',NC,' 'JNK,'] Model Output']) 
grid 
Figure(15) 
plot(difft,'r’) 
axis({5000,5 100,-0.5,0.5]) 
title({'Actual Output - ARMAX[',NA,'’,NB,' ',NC,' '",NK,']Model Output *Data 
Filtered*']) 


grid 
avediff=mean(abs(diff)) % Ave. output difference 
avedft=mean(abs(difft)) % Ave. output difference (Prefiltered) 


% Generate Zero-Pole Plots 

Figure(12) 

zpplot(th2zp(arxm),3,[],1.2) 

title({"Zero Pole Plot ARMAX MODEL|',NA,’ 'JNB,' '"JNC,' ',NK,']']) 

Figure(13) 

zpplot(th2zp(arxmft),3,[],1.2) 

title(("Zero Pole Plot ARMAX [',NA,'',NB,’ ',NC,''",NK,'] Model *Data Prefiltered*']) 


To To To To Yo Yo To Yo To To Yo Yo To To Yo To To To Yo To Fo To To To To To To To Yo To Yo Fo Yo To Yo To Yo To Yo To Vo Vo Vo 


% bjmod.m %o 
% Box-Jenkins Model Construction and Analysis Jo 
%o George D. Beavers Jo 
% June 1997 Wy 


Jo To To To To To To To To To To To Vo To To To To To To To To To To Mo To Yo To Mo To To To To To To To Yo To To To Yo To Yo Lo 
% This program generates a model based on the input-output data of asystem using % 
% the Box-Jenkins parameter estimation method. It also generated a series of graphical% 
% and numerical results in order to evaluate the suitability of the model %o 
Jo To To To To To To To To To To To To To To To To Yo To Mo To To Vo To To Yo Vo Vo Yo To To To Yo Yo To Yo To Vo To Yo To To Fo 
function [bjm,bjmft]=bjmod(j,nns,itr,z2,z3) 

To To To To To Yo To To Yo To To To To Mo To To Yo Yo To To Vo Yo Yo Vo Yo Vo 


ZY bjym=Box-Jenkins based model %o 
% bjmft=Box-Jenkins based model (filtered data) % 
% j=actuator %o 
Yo nns=structure([na nb nc]) % 
% itr=number of iterations To 


20 


% z2=input-output data for model construction % 
% Z3=input-ouput data for model evaluation % 
JoJo To To To Lo To Yo Yo Yo Lo Yo To To To Yo Yo To Yo To Yo Yo Yo To Yo Vo 
nb=nns(1,1); 

NB=int2str(nb); 

nc=nns(1,2); 

NC=int2str(nc); 

nd=nns(1,3); 

ND=int2str(nd); 

nf=nns(1,4); 

NF=int2str(nf); 

nk=nns(1,5); 

NKe=int2str(nk); 


% This section loads separate data to compute the actual actuator frequency response % 


% and then plot it against the model’s frequency response % 
arxm=bj(z2,nns, itr); % Build the model % 
arxm=sett(arxm,0. le-3); % Set the sampling interval % 
garxm=th2ff(arxm); % Convert to frequency function % 

eval([‘load s',int2str(q),'m']) % Load data for actuator frequency function% 
X1=0211x; 

Y 1=20.*log10(abs(0211)); 

Figure(2) 

semilogx(X1,Y1,'b’), 

hold on 


semilogx(garxm(:,1)/(2*p1),20*log lO(garxm(:,2)),'r'); 
title(["Blue=Actuator#',num2str(j),, Red=Box-Jenkins MODEL|',NB,' ',NC,' ', ND,’ ',NF, 
SNK, '}']) 

xlabel(‘Frequency (Hz)’) 

ylabel(‘Magnitude(dB)’) 

grid 

hold off 


eval([‘load s',int2str(j),'p'J) 

X2=0 211%: 

Y2=180*unwrap(angle(o211))/pi; 

Figure(3) 

semilogx(X2, Y2,'b') 

hold on 

semilogx(garxm(:,1)/(2*p1),(garxm(:,4)),'r'); 

title({'Blue=Actuator#' num2str(j),,Red=Box-Jenkins MODEL[',NB,' ',NC,'',ND,' '",NF, 
ENED 


xlabel(‘Frequency (Hz)’) 
ylabel(‘Phase (Deg)’) 
grid 

hold off 


J The above process is repeated using prefiltered data. A 10" order Butterworth filter % 
% is used. The passband is given in terms of fractions of the Nyquist frequency (0.52)% 


% This is a low pass filter with an upper bound of approximately 1|.5kHz J 
Ze—romlt(z2,10,.52); % Prefilter model construction data % 
arxmft=bj(zf,nns, itr); % Construct model % 
arxmf{t=sett(arxmft,0. le-3); % Set sampling interval % 
garxmft=th2ff(arxmft); Y Convert to frequency function % 

eval([‘load s',int2str(j), m']) % Load data for actuator frequency response 
Pol=0211Xx: 


Y 1=20.*log 10(abs(0211)); 


Figure(4) 

semilogx(X1,Y1,'b’) 

hold on 

semilogx(garxmft(:, 1)/(2*p1),20*log 10(garxmft(:,2)),'r'); 

title({‘Blue=Actuator#' num?2str(j),,Red=Box-Jenkins MODEL[’,NB,' ',NC,' ',ND,' ',NF,’ 
',NK,']*Data Prefiltered*']) 

xlabel(‘Frequency (Hz)’) 

ylabel(‘Magnitude(dB)’) 

grid 

hold off 


eval([‘load s',int2str(j),'p']) 
x2 =0211x: 
Y2=180*unwrap(angle(o211))/pi; 


Figure(5) 

semilogx(X2, Y2,'m’) 

hold on 

semilogx(garxmft(:, 1)/(2*p1),(garxmft(:,4)),'g'); 

title((‘Blue=Actuator#' ,num2str(j),,Red=Box-Jenkins MODEL[',NB,'';NC,’ ',ND,' ',NF,. 
',NK,']*Data Prefiltered*']) 

xlabel(’‘Frequency (Hz)’) 

ylabel(’Phase (Deg)’) 

grid 

hold off 


10] 


% ‘This section generates plots and numerical results used in evaluating the model % 
Y Comparison of model output to actual output % 
% Apply the input row of the validation data to the model and then compare with Jo 
% the actual output of the system corresponding to the validation data and plot results 

J 


Figure(6) 

{yh,fit]=compare(z3,arxm); 

axis({5000,5 100,-1.0,1.0]) 

title({'Plot of Box-Jenkins [',NB,'',NC,' ',ND,' ',NF,'',NK,'] Output vs. Actual output’ ]) 
text(5020,0.75,'Red:Model Output, Blue: Measured Output’) 

xlabel(‘Sample’) 


Figure(7) 

[yhf,fitf]=compare(z3,arxmft); 

title({'Plot of Box-Jenkins MODEL [',NB,'',NC,'',ND,'',NF,' ',NK,'] Output vs. Actual 
output *Data Prefiltered*')) 

text(5020,0.75,'Red:Model Output, Blue: Measured Output’) 
axis({5000,5100,-1.0,1.0]) 

xlabel(‘Sample’) 


% Calculate residuals and auto and cross correlation functions and plot results 


Figure(8) 

el=resid1(zfv,arxm); 

Figure(9) 

e2=resid1(zfv,arxmft); 

Figure(10) 

plot(el,'m'),grid 

title({"Residual of Box-Jenkins Model [',NB,'',NC,' ',ND,' ',NF,' ',NK,'J']) 
ylabel(‘Residual’) 

xlabel(‘Sample’) 

axis({5000,5 100,-0.1,0.1]) 

Figure(1 1) 

plot(e2,'m’'), grid 

title({"Residual of Box-Jenkins Model [',NB,'',NC,' '",ND,' ',NF,'',NK,'] *Data 
Prefiltered*']) 

ax1s({[5000,5 100,-0.05,0.05]) 

xlabel(‘Sample’) 

ylabel(’Residual’) 
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% Numerical results 
el pos=abs(e1); 
e2pos=abs(e2); 
ypos=abs(z3(:,1)); 


ave=mean(e 1 pos) % Average residual 
aveft=mean(e2pos) % Average residual (Prefiltered Data) 
outave=mean(ypos) % Average actuator output 

fit Yo Mean Square Fit calculated above 
fitf % Mean Square Fit (Prefiltered) 


% Calculate and plot actual output minus model output 
Zooul—75(:,1); 
for i=5000:5100; 
difft(i)=z3out(i)-yhf(i); 
end 
Zoout=Z5(:,1); 
for j=5000:5 100; 
diff(j)=z3out(j)-yhG); 
end 
Figure(12) 
plot(diff,'r') 
ax1S([5000,5 100,-0.5,0.5}) 
title(‘Actual Output - Model Output’) 
Figure(13) 
plot(difft,'r’) 
ax1s({[5000,5 100,-0.5,0.5]) 
title(‘Actual Output - Model Output *Data Filtered*') 
avediff=mean(abs(diff)) % Average output difference 
avedft=mean(abs(difft)) % Average output difference(Prefiltered) 


% Generate Zero-Pole Plot 

Figure(14) 

zpplotl(th2zp(arxm),3,[],1.2) 

title(['Zero Pole Plot Box-Jenkins Model [',NB,'',NC,'',ND,' '|NF,' '|NK,'} ']) 
Figure(15) 

zpplotl(th2zp(arxmft),3,[], 1.2) 

title([‘Zero Pole Plot Box-Jenkins Model [',NB,'',NC,'',ND,' '"JNF,'',NK,'] *Data 
Eretilicred= |) 


Jo To To To Fo Po To Yo To Po To To To Yo Yo Fo To To To To To To Fo Fo To To Fo To To To To To To Yo To Vo Vo Fo To To To Yo Yo 


% plotarx.m Jo 
Jo Compare Three ARX Model Structures Jo 
To George D. Beavers Jo 
Jo iiittemoo7 Jo 


Jo To To To To To To To To Wo Yo Fo Yo Vo Vo To Yo To To Vo To Lo To To Fo Po Fo To To To To To To To Fo To Fo To To To To Fo To 
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% This program generates the frequency response plots of the three best ARX model % 
% structures for evaluation and structure selection Jo 
ToT To To To To To To To Yo To To To Yo Yo Yo To To To To To To To Yo To To To Yo To To To To To Yo To Vo To Yo To To Yo Yo Yo 


j=input(‘enter strut>>'); 


eval([‘load 4h',int2str()]) 
eval([‘load s’,int2str(),'m']) 


z2=[trace_y(j,1:10000)' trace_y(7,1:10000)']; 
z2=dtrend(z2); 

z3=[trace_y(j, 10001:20000)' trace_y(7,10001:20000)'}; 
z3=dtrend(z3); 


nnl=[7 6 1] % Top three structures 
nnZ=> 1S) 
iiss levee 


X1=02i1x; _ 
Y 1=20.*log10(abs(0211)); 


zf=idfilt(z2, 10,.52); 

% Construct the models % 
fml=arx(zf,nn 1); 
fm1=sett(fm1,0.1le-3); 
efml=th2ff(fm1); 


fm2=arx(zf,nn2); 
fm2=sett(fm2,0. le-3); 
gfm2=th2ff(fm2); 


fm3=arx(zf,nn3); 
fm3=sett(fm3,0. le-3); 
efm3=th2ff(fm3); 


% Generate comparison plots 
Figure(1) 

semilogx(X1,Y1,'b'), 

title([‘Strut #',int2strq),’ Magnitude’]) 
xlabel(‘Frequency (Hz)’) 
ylabel(Magnitude(dB)’) 

hold on 
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semilogx(gfm1(:,1)/(2*pi),20*log10(gfm1 (:,2)),'r'); 
semilogx(gfm2(:,1)/(2*pi),20*log 10(gfm2(:,2)),'m'); 
semilogx(gfm3(:,1)/(2*pi),20*log10(gfm3(:,2)),'w');,grid; 
hold off 


Fioure(2); 

eval([‘load s',int2str(j),'p']) 

X2=0211x; 
Y2=180*unwrap(angle(o2i1))/pi; 
semilogx(X2,Y2,'b'), 

title([‘Strut #',int2strG),’ Phase']) 
xlabel(‘Frequency (Hz)’); 

ylabel(‘Deg’); 

hold on 
semilogx(gfm1(:,1)/(2*pi),(gfm1(:,4)),'r'); 
semilogx(gfm2(:,1)/(2*pi),(gfm2(:,4)),'m’); 
semilogx(gfm3(:,1)/(2*p1),(gfm3(:,4)),'w');,grid; 


hold off 

To To To To To To To To To To To Yo Fo To To To To Yo Yo Fo To Yo Lo To To To To Mo To To To To Yo To To Yo To To To To Yo Yo To 
% plotamb.m % 
% Model Comparison % 
% George D. Beavers % 
Jo June 1997 Jo 


To To To To To To To To To To To To To To To To To To To To To To To Yo Yo To To Yo To To Yo Yo To To To To To To To Yo Yo To Yo 
% This program generates the frequency response plots of the three best models based % 
7 on the ARX structure[7 6 1]. ARX, ARMAX and BOX-JENKINS models are % 
% compared to select the best model Jo 
Jo To To To To To To Yo To To To To To To To To To To To To To Zo To To To To Fo To Yo Yo To To To To To To To To To To Yo Yo Yo 
j=input(‘enter strut>>’); 


eval({‘load 4h’,int2str(j)]) 
eval([‘load s’,int2str(j),'m']) 


z2=(trace_y(j,1:10000)' trace_y(7,1:10000)']; 
z2=dtrend(z2); 

z3=[trace_y(j,10001:20000)' trace_y(7,10001:20000)']; 
z3=dtrend(z3); 


% Define model structure % 
fiie—=(7 62 |] 

nnz—(/ 6 1} 

Nie —|Oree ay all 

Dl=o21lx: 
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Y1=20.*log10(abs(02i1)); 
% Generate models (non-filtered) % 
%m1=armax(z2,nn 1,60); 
%m1=sett(m1,0.1le-3); 
%gm1=th2ff(m1); 
%m2=arx(Z2,nn2); 
%m2=sett(m2,0.1le-3); 
%ogm2=th2ff(m2); 
%om3=bj(z2,nn3,60); 
%Jm3=sett(m3,0.1le-3); 
%gm3=th2ff(m3); 


zf=idfilt(z2,10,.52); 


% Generate Models (Prefiltered) 
fml=armax(zf,nn1 ,60); 
fm1=sett(fm1,0.1le-3); 
gfm1=th2ff(fm1); 


fm2=arx(zf,nn2); 
fm2=sett(fm2,0.1le-3); 
gfm2=th2ff(fm2); 


fm3=b)(zf,nn3,60); 
fm3=sett(fm3,0. 1le-3); 
gfm3=th2ff(fm3); 


% Generate plots 

Figure(1) 

semilogx(X1, Y1,'b') 

title([‘Strut #',int2strq),, Magnitude’}) 
xlabel('Frequency (Hz)’) 

ylabel('Magnitude(dB)’) 

hold on 

semilogx(gfm1 (:,1)/(2*pi),20*log10(gfm1(:,2)),'r'); 
semilogx(gfm2(:,1)/(2*pi),20*log 10(gfm2(:,2)),'m’); 
semilogx(gfm3(:,1)/(2*pi),20*log 10( gfm3(:,2)),'w');,grid; 
hold off 


Figure(2); 

eval([‘load s',int2str(j),'p']) 
AZ=OlT IE 
Y2=180*unwrap(angle(o211 ))/p1; 
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semilogx(X2, Y2,'b') 

title([{‘Strut #',int2str(),’ Phase')) 
xlabel('Frequency (Hz)’); 

ylabel('‘Deg’); 

hold on 
semilogx(gfm1(:,1)/(2*p1),(gfm 1(:,4)),'r'); 
semilogx(gfm2(:,1)/(2*pi),(gfm2(:,4)),'m’); 
semilogx(gfm3(:,1)/(2*p1),(gfm3(:,4)),'w’);,grid; 


hold off 

Jo To To To To To To Lo To To To To Fo To Yo To To Yo Lo To Lo To To Fo To To To Yo To Yo To Yo No Yo Fo To To Yo To Yo To No Yo 
% compid.m %o 
Jo Interaction Comparison Yo 
Jo June 1997 Jo 


To To To To To To Lo Fo To Lo To Yo Yo Lo Yo To To To To To To Lo To Lo To To To To Lo Yo To Yo Vo To To Lo Yo Lo To To Yo Vo Yo 
% This program generates the frequency response plots of the actuators when actuator % 
% number one is excited with 50mV “pink noise”. These plots use the ARX[7 6 1] Jo 


% model frequency repines of the respective actuator. This is compared with the Jo 
% modeled response of actuator number one These plots are used to check for Jo 
% interaction between actuators Jo 


To To To To To To Lo To To To To Yo Yo Lo To To To Lo Lo To Lo Lo To Fo To To Lo To To To Fo To To Fo Lo Lo Yo To Fo To Yo Vo Yo 


function []=compid(j) 
% j= actuator number% 
int2str(j) 


disp([‘load 4h’ int2str(j)}) 
eval([‘load ' '4h' int2str(j) '.mat']) 


z1=[trace_y(1,1:10000)' trace_y(7,1:10000)'}; 
Zl—atrematz ||); 

zla=[trace_y(1,10001:20000)' trace_y(7,10001:20000)'); 
Zla=dtrend(zla); 


ml—=arx(Z ie) ot}: 


th1=sett(th 1,0. le-3); 
eth |=th2ff(th 1); 


[OF 1=2:6 
Z2=[trace_y(i,1:10000)' trace_y(7,1:10000)']; 


Zo=dtrenal ZZ): 
z7=(trace_y(i,10001:20000)' trace_y(7,10001:20000)'}; 
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z/=dtrend(z7/); 


th=arx(Z2.17 GH): 
th=sett(th,0.1e-3); 
eth=th2ff(th); 


Figure(i-1); 


semilogx(gth1(:,1)/(2*p1),20*log 10(gth1(:,2)),'b'); 

hold on 

semilogx(gth(:,1)/(2*pi),20*log10(gth(:,2)),'r'); 

title Magnitude: Input=50mV Pink Noise’) 

xlabel('Frequency (Hz)’) 

ylabel('dB') 

text(10,45, [Excite Actuator # num2str(j) ‘Blue; Sense Actuator ', 
num2str(i),=Red']) 


erid 
hold off 
end 
for 1=2:6 


z2=([trace_y(i, 1:10000)' trace_y(7,1:10000)'); 
z2=dtrend(z2); 

z7=[trace_y(i, 10001:20000)' trace_y(7,10001:20000)'); 
z/7=dtrend(z7); 


th=arx(z2,[7 6 1]); 
th=sett(th,0.1e-3); 
eth=th2ff(th); 


Figure(i+4) 

semilogx(gth1(:,1)/(2*pi),(gth1(:,4)),'b’); 

hold on 

semilogx(gth(:,1)/(2* p1),(gth(:,4)), 1’); 

title(/Phase: Input=5Omv Pink Noise’) 

text(10,-550,[Excite Actuator #' num2str(j) '=Blue; Sense Actuator ',snum2str(i), 
=Red}) 

xlabel(’‘Frequency (Hz)’) 

ylabel(‘Degrees') 

grid 

hold off 
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end 


Jo To To To To To To To To To Yo To To Fo Yo To Fo To Lo To To Yo To To To To Yo To To Yo Yo Yo Yo To Yo Mo To Yo Yo To Yo To Zo 


% fig1.m, fig2.m, fig3.m, fig4.m Yo 
%o PSD of Applied Controllers % 
% June 1997 %o 
ToT To To To To To Fo 0 To Yo Yo To To Vo Fo To Vo Lo To To Yo To To To To To To To Yo Yo To Yo Yo To Yo To Yo Yo Vo Yo Yo Yo 
% This program generates the power spectral density plots resulting from applying % 
% a lead-lag compensator to an individual actuator and the all six actuators to evaluate % 
% the effects of actuator interaction Jo 
To To To To To To To Vo Yo To To Yo Lo To Yo Yo To Yo To To Yo To To To Yo To Yo Yo Vo To Fo To To Yo To Yo Yo To To To Yo Yo Yo 
%fig1.m7% 

% Plots over range of 0-5kHz% 

clear 

load s1 Yo Open loop data with 42 Hz disturbance on% 

Figure(1) 

clg 


psd(trace_y(2,1:20000), 1024*1,10000,[]); 

%[popen,fopen]=psd(trace_y(2,1:20000), 1024*12,10000,[]); 

hold on 

load s2 %oCompensator applied to actuator #2 ( 

psd(trace_y(2, 1:20000),1024*1,10000,kaiser(512,5),256); 
%o[ps2,fs2]=psd(trace_y(2,1:20000), 1024*12,10000,[]); 

load s3 %Compensator applied to all six actuators 

psd(trace_y(2, 1:20000), 1024*1,10000,[]); 

%o(ps3,fs3 ]=psd(trace_y(2,1:20000), 1024*12,10000,[]); 

hold off 
ZYosemilogx((fopenk),20*log 1 O(popenk),'r',(fs2k),20* log 10(fs2k),'y',(fs3k),20*log 10(fs3k) 
0) 

Title(‘Blue-Open Loop; Magenta-Strut #2 Closed Loop; Red-All 6 Struts Closed Loop’ ) 


To To To To To To To To To To To To To To Fo To To Vo To To To To To To To To To To Yo To To To Uo To To Vo Fo Fo To Yo Yo Yo Vo 
%fig2.m% 
% Plots same Figure as fig1.m over range of 0-200 Hz% 


clear 

load s1 %Open loop-42 Hz Disturbance on% 
Figure(2) 

clg 

psdb(trace_y(2,1:20000),1024*8, 10000,[]); 

%o| popen,fopen]=psd(trace_y(2, 1:20000),1024* 12, 10000,[]); 

hold on 
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load s2 %Compensator applied to actuator #2% 
psd(trace_y(2,1:20000), 1024*8,10000,[]); 

%(ps2,fs2]=psd(trace_y(2,1:20000), 1024*12,10000,[]); 

load s3 Y% Compensator applied to actuator #2% 
psd(trace_y(2,1:20000), 1024*8, 10000,[]); 

%[ps3,fs3 ]=psd(trace_y(2,1:20000), 1024* 12,10000,[]); 

hold off 

%osemilogx((fopenk),20*log 10(popenk),'r'’,(fs2k),20*log 10(fs2k),'y',(fs3k),20*log 10(fs3k) 
’) 

Title(Blue-Open Loop; Magenta-Strut #2 Closed Loop; Red-All 6 Struts Closed Loop’ ) 
% Title(‘(Open Loop-42Hz Disturbance’ ) 

axis({O 200 -70 30]) 


To To To To To To To To To To To Yo To To Yo Yo Go Yo Yo Yo To Yo Yo Yo Yo Yo To To To Yo Lo No To To To To To To To To Yo To Vo 
%ofig3.m% 
Yplots PSD sensed by 3 geophones; actuator #2, #3 and #6 over range of 0-200Hz% 


clear 

load s2 % Compensator applied to actuator #2 % 
Figure(3) 

clg 

psd(trace_y(2,1:20000), 1024*1,10000,[]); 
%(popen,fopen]=psd(trace_y(2, 1:20000), 1024*12,10000,[]); 

hold on 

psd(trace_y(3,1:20000), 1024*1,10000,[]); 
%[ps2,fs2]=psd(trace_y(3,1:20000), 1024*12,10000,[]); 
psd(trace_y(6,1:20000), 1024* 1,10000,[]); 

% [ps3 ,fs3]=psd(trace_y(6,1:20000), 1024*12,10000,[}); 

hold off 

%semilogx((fopenk),20*log 1O(popenk),'r',(fs2k),20*log 1 0(fs2k),'y',(fs3k),20*log 1 0(fs3k) 
‘) 

Title(‘Blue-#2 Geophone; Magenta-#3 Geophone; Red-# 6 Geophone' ) 
%axis({O 200 -70 30]) 


Zo To To To To To To To To To Yo To Yo Yo To To Yo To Yo Yo To To To Yo Yo To Lo To To Yo Yo To To To To Yo To Yo To To To To Yo 
%fig3.m% 
Yoplots PSD sensed by 3 geophones; actuator #2, #3 and #6 over range of 0-5kKHz% 


clear 

load s2 %oCompensator applied to actuator #2 % 
Figure(4) 

clg 

psd(trace_y(2,1:20000), 1024*8,10000,[]); 
%o{popen,fopen]=psd(trace_y(2, 1:20000), 1024* 12, 10000,[]); 

hold on 
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psd(trace_y(3,1:20000),1024*8,10000,[]); 

Jo[ps2,fs2]=psd(trace_y(3,1:20000), 1024*12,10000,[]); 

psd(trace_y(6, 1:20000), 1024*8, 10000,[]); 

Jo ps3,fs3]=psd(trace_y(6,1:20000), 1024* 12, 10000,[]); 

hold off 

Jsemilogx((fopenk),20*log10(popenk),'r' ,(fs2k),20*log10(fs2k),'y',(fs3k),20*log 1 O(fs3k) 
b) 

Title(Blue-#2 Geophone; Magenta-#3 Geophone; Red-# 6 Geophone' ) 

axis({O 200 -70 30]) 
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APPENDIX B. MODEL STRUCUTURE SELECTION PLOTS 
A. Structure Acceptance Plots for Actuator Two 


The following plots were utilized in selecting the ARX structure [na=7 nb=6 nk=1] 
for actuator number two. 
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Frequency Response Comparison 
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Figure A.7 
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B. Model Structure Acceptance for Actuator Three 


The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=1] 
for actuator number three. 
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C. Model Structure Acceptance for Actuator Four 


The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=1] 
for actuator number four. 
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D. Model Structure Acceptance for Actuator Five 


The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=1] 
for actuator number five. 
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E. Model Structure Acceptance for Actuator Six 


The following plots were utilized in selecting the ARX structure [na=7 nb=7 nk=1] 
for actuator number six. 
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APPENDIX C. MODEL VALIDATION PLOTS 
A. Application of ARMAX Method to Actuator Number One. 
The following plots were used in evaluating the application of the ARMAX and 


Box-Jenkins parameter estimation methods in modeling actuator number one. The 
structure [na=7 nb=6 nc=2 nk=1] was used. 
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Frequency Response Comparison 
Plot of ARMAX MODEL? 6 2 1] Output vs. Actual output *Data Prefiltered* 
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Figure A.6 
Zero-Pole Plot 
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B. Application of ARMAX Method to Actuator Number Two. 


The following plots were used tn evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number two. The structure [na=7 nb=6 
nc=2 nk=1] was used. Two additional plots from the Box-Jenkins [nb=6 nc=2 nc=2 na=7 
nk=1} model are tncluded for comparison. 
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Figure B.7 
Zero-Pole Plot for Box-Jenkins Model 
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C. Application of ARMAX Method to Actuator Number Three. 


The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number three. The ARX based 
structure [na=7 nb=6 nc=2 nk=1] was used. 
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Figure C.6 
Zero-Pole Plot 
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D. Application of ARMAX Method to Actuator Number Four. 


The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number four. The structure [na=7 nb=6 


nc=2 nk=1] was used. 
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Frequency Response Comparison 
Plot of ARMAX MODEL? 6 2 1] Output vs. Actual output 
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EK. Application of ARMAX Method to Actuator Number Five. 


The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number five. The structure [na=7 nb=6 
nc=2 nk=1] was used. 
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F. Application of ARMAX Method to Actuator Number Six. 


The following plots were used in evaluating the application of the ARMAX 
parameter estimation method in modeling actuator number six. The structure[na=7 nb=6 


nc=2 nk=1] was used. 
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Frequency Response Comparison 
Plot of ARMAX MODEL[7 6 2 1] Output vs. Actual output *Data Prefiltered* 
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